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ABSTRACT 


This  thesis  presents  a  computationally-oriented 
linear  boundary  collocation  technique  for  obtaining  a 
close  approximation  to  the  solution  of  the  Dirichlet 
problem  for  Laplace's  equation.  Necessary  concepts  and 
algorithms  from  the  theory  of  linear  approximation  are 
reviewed,  and  numerical  results  are  given. 
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CHAPTER  I 


INTRODUCTION 


This  thesis  reviews  the  exchange  algorithm  for  computing 
a  best  linear  approximation  on  a  discrete  point  set  and  shows 
how  it  may  be  used  by  an  algorithm  of  Remes  to  compute  a 
best  linear  approximation  on  a  continmum.  These  two 
algorithms  are  employed  by  a  linear  boundary  collocation 
method  to  compute  a  close  approximation  to  the  solution  of 
the  Dirichlet  problem  for  Laplace's  equation. 

1.1  The  Dirichlet  Problem  for  Laplace's  Equation 

Consider  the  following  classical  problem  from  the 
theory  of  partial  differential  equations:  Let  C  be  a 
simple  closed  curve,  R  the  region  interior  to  C  ,  and 
f  a  continuous  function  on  C  .  Find  a  function  u  *  u(x,y) 
having  continuous  second  order  derivatives  in  R  satisfying 


(1.1) 


and 


(1.2) 


u  =  f  on  C 


The  subscripts  in  (1.1)  denote  partial  differentiation. 
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Equation  (1.1)  is  known  as  Laplace’s  equation.  Functions 
having  continuous  second  order  derivatives  in  R  and 
satisfying  (1.1)  are  called  harmonic  functions  in  R  . 

Equation  (1.2)  is  the  Dirichlet  boundary  condition  for 
Laplace’s  equation.  Determining  the  solution  u  is  the 
Dirichlet  problem  for  Laplace's  equation.  Ahlfors 
([1],  PP*  175-199)  has  proven  that  u  exists  and  is  unique. 

The  Dirichlet  problem  arises  in  a  wide  variety  of 
physical  contexts,  including  heat  transfer,  hydrodynamics, 
electricity,  and  aerodynamics.  Sample  problems  may  be 
found  in  Forsythe  and  Wasow  [16]  and  Fox  [17]. 

Except  for  some  simple  boundaries  such  as  the  rectangle 
and  the  circle,  where  analytical  solutions  can  be  obtained 
in  closed  form,  most  boundaries  require  that  some  numerical 
method  be  used  to  obtain  the  solution.  Most  of  the 
numerical  work  devoted  to  the  Dirichlet  problem  has  been 
concerned  with  replacing  the  derivatives  in  (1.1)  by 
appropriate  finite  differences  and  then  determining  the 
solution  at  a  finite  set  of  points  in  the  region  R  .  How¬ 
ever,  this  technique  is  difficult  to  apply  to  regions  whose 
boundaries  are  irregular,  and  the  resulting  table  of 
solution  values  is  not  convenient  for  use  in  subsequent 
numerical  work.  Attempts  have  therefore  been  made  to  derive 
methods  which  will  not  only  yield  the  solution  in  closed 
form,  but  also  be  easy  to  apply  to  irregularly  shaped  boundaries. 
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1 . 2  Collocation  Methods 

Attempts  to  solve  boundary-value  problems  in  partial 
differential  equations  are  often  made  using  collocation 
methods.  One  assumes  as  an  approximation  to  the  solution 
u(x,y)  an  expression  of  the  form  g(x,y;c)  which  depends 
on  a  vector  of  parameters  c  =  (c^#...,c  )  such  that  for 
arbitrary  values  of  c 

a)  the  differential  equation  is  already  satisfied 
exactly,  or 

b)  the  boundary  conditions  are  already  satisfied  exactly, 
or 

c)  g  satisfies  neither  the  differential  equation  nor 
the  boundary  conditions. 

One  then  tries  to  determine  c  so  that  g  satisfies 
in  method  (a),  the  boundary  conditions,  or 
in  method  (b),  the  differential  equation,  or 
in  method  (c),  the  boundary  conditions  and  the  differen¬ 
tial  equation 

as  accurately  as  possible  in  some  sense  yet  to  be  defined. 

Definition  1.1.  A  collocation  method  of  type  (a) 
described  previously  is  called  a  boundary  collocation  method. 

For  partial  differential  equations,  boundary  collocation 
methods  are  usually  preferred  for  two  reasons.  First,  it  is 
generally  easier  to  find  functions  which  satisfy  a  particular 
differential  equation  than  it  is  to  find  functions  which 
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satisfy  the  boundary  conditions  for  a  particular  boundary. 
Second,  boundary  collocation  methods  are  usually  chosen 
since  their  use,  in  so  far  as  integration  is  concerned, 
requires  the  evaluation  of  integrals  around  the  boundary  C 
rather  than  throughout  the  region  R  .  These  observations 
were  also  verified  in  practice  by  Shuleshko  [29],  who  made 
a  comparative  analysis  of  the  three  collocation  techniques 
applied  to  the  torsion  problem  and  concluded  that  not  only 
were  boundary  collocation  methods  easier  to  apply,  but  also 
yielded  better  results. 

Linear  partial  differential  equations  have  the  property 
that  au  +  3v  is  a  solution  of  the  partial  differential 
equation  if  u  and  v  satisfy  the  equation  and  a  and  3 
are  any  two  constants.  Hence  boundary  collocation  methods 
for  linear  partial  differential  equations  can  have  a 
special  form. 


Definition  1.2.  A  boundary  collocation  method 
for  a  linear  partial  differential  equation  which  uses  an 


n 


approximating  function  of  the  form  g  =  £  c  .g  .  ,  where 

tJ  J 


Q-l 


each  function  g.  ,  j=l3...3n  ,  satisfies  that  partial 
differential  equation ,  is  called  a  linear  boundary 


collocation  method. 
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1.3  Linear  Boundary  Collocation  Methods  for  the  Dirichlet 
Problem 


Because  Laplace’s  equation  is  linear,  and  since  the 

sum  of  two  continuous  functions  is  continuous,  the  following 

is  a  linear  boundary  collocation  method  for  the  Dirichlet 

problem:  Let  g^»t..#g  be  harmonic  functions  in  the 

region  R  .  Then  for  an  arbitrary  vector  of  coefficients 

n 

c  =  (c^,...,cn),  J  cjsj  is  also  a  harmonic  function  in 
R  .  All  linear  boundary  collocation  methods  are  dis¬ 
tinguished  by  the  manner  in  which  c  is  determined. 


1.3.1  Interpolation  Techniques  In  this  method,  n 

distinct  points  s^,...,s  are  selected  on  the  boundary 

n 

C  ,  and  the  parameters  are  chosen  so  that  £  c . g . 

j=l  J  J 

interpolates  to  f  at  these  points.  That  is,  (c^,...,c  ) 
is  the  solution  of  the  system  of  linear  equations. 


n 

(1.3)  l  Cjgj(si)  =  f ( s± )  ,  i=l , . . . , n  . 

J  —  b 


The  essential  idea  of  this  technique  was  first  suggested 
by  Walsh  [34]  in  1929,  when  he  discussed  the  problem  of 
approximating  harmonic  functions  by  harmonic  polynomials. 
Since  then  the  method  has  been  applied  with  varying 
success  to  a  number  of  engineering  problems,  including  ones 
in  heat  transfer,  elasticity,  and  the  bending  and  buckling 
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of  plates  [^,5,6,7,8,19,23,26,28,30,33].  The  technique 
was  first  dubbed  the  "point-matching  method"  by  Conway  [6] 
in  1961. 

Not  too  much  is  known  about  this  method  theoretically. 
Its  use  seems  to  indicate  that  its  success  depends  upon  the 
geometry  of  the  problem,  but  this  dependence  has  not  been 
explicitly  determined.  Another  problem  arises  from  the 
fact  that  the  coefficient  matrix  in  (1.3)  may  be  singular 
for  some  sets  of  points  on  C  .  Further,  it  is  not  known 
under  what  conditions  the  method  converges.  Perhaps  the 
most  critical  problem  is  that  a  general  set  of  properties 
of  the  "best"  set  of  matching  points  has  not  been  found. 

Theoretical  material  concerning  this  technique  is 
discussed  in  a  series  of  papers  by  Curtiss  [9-11]  and 
Sobczyk  [31].  On  the  basis  of  simple  cases  of  circles  and 
ellipses  that  have  been  investigated  successfully  [35,36], 
and  also  arguing  by  analogy  from  the  case  of  complex 
analytic  interpolation,  Curtiss  [9]  concludes  that  for 
sufficiently  regular  boundaries  there  probably  exist 
sequences  of  points  for  which  this  interpolation  process 
converges  for  a  wide  class  of  boundary  data.  However,  his 
suggestion  that  these  points  are  probably  obtainable  as 
the  exterior  conformal  image  of  equally  spaced  points  on 
the  unit  circle  (such  as  the  roots  of  unity)  is  unsuitable 
for  numerical  work  because  the  conformal  map  is  not  always 


. 


’ 

•  L  .  • 


7 


available.  Sobczyk  [31]  has  even  shown  that  the  resulting 
image  points  on  C  may  result  in  a  singular  coefficient 
matrix  in  (1.3). 

1.3*2  Least-Squares  Techniques  In  a  least-squares 
approach,  c  is  chosen  to  minimize 

(1.4)  E,  ■  /  [f (s)  -  l  c.g,(s)]2  ds  . 

1  C  j=l  J  J 

Usually  the  integration  in  (1.4)  cannot  be  carried  out 
analytically.  Therefore,  m  points  are  chosen  on  the 
boundary  C  and  m  weights  w^  corresponding  to  a 
quadrature  rule,  and  then  the  vector  c  is  chosen  to 
minimize 


(1.5) 


m  n 

E2  =  w  1[f(s1)  -  ^  cjgj(8i)] 


In  contrast  to  the  interpolation  technique,  the  discrete 
process  in  (1.5)  should  be  convergent  as  n  becomes  large 
and  m>>n  .  Details  of  this  method  and  numerical  results 
can  be  found  in  papers  by  Lo  [20],  Davis  [12],  Davis  and 
Rabinowitz  [14],  O'Jalvo  and  Linzer  [24],  Merriman  [21], 
and  Sparrow  [32]. 


. 
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1 . 4  Purpose  of  Study 

In  Chapter  II,  we  shall  review  some  basic  theory 
relevant  to  the  Dirichlet  problem  and  thus  attempt  to  show 
that  choosing  c  to  minimize 

n 

(1.6)  M  =  max  |f(s)  -  J  c.g(s)| 

seC  j  =  1  J  J 

is  a  more  natural  linear  boundary  collocation  method  than 
either  the  point-matching  method  or  the  least-squares 
method.  In  Chapter  III,  IV  and  V  we  shall  show  that  such 
a  vector  c  always  exists  and  we  shall  present  a 
computationally  oriented  algorithm  for  finding  it.  Results 
of  numerical  computations  will  be  given  in  Chapter  VI,  and 
conclusions  and  suggestions  for  further  research  in 
Chapter  VII. 

A  discrete  version  of  this  approach  has  been  made  by 
Quon  and  Leung  [27]  using  linear  programming  techniques. 


. 
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CHAPTER  II 


BASIC  THEORY 


2 . 1  Some  Concepts  in  Linear  Approximation  Theory 

The  problem  outlined  in  the  previous  chapter  (see 
Section  1.4)  is  a  typical  one  in  approximation  theory  in 
that  it  involves  the  selection  from  a  given  set  of  functions 
one  that  is  in  some  sense  "close"  to  a  prescribed  function 
(usually)  not  in  the  set.  Any  discussion  of  approximation 
requires  that  we  have  a  way  of  measuring  the  "closeness" 
of  two  functions.  It  is  reasonable  to  demand  that  this 
measure  should  have  some  of  the  properties  of  geometrical 
distance,  or  more  precisely,  that  it  should  be  a  metric. 

2.1.1  Metric  Spaces 

Definition  2.1.  Let  S  be  a  set  of  elements 
and  let  d  be  a  real-valued  function  defined 
on  pairs  of  elements  of  S  .  Then  d  is  called  a  metric 
if  it  possesses  the  following  properties: 

i)  d(xsx)  =  0 , 
ii)  d(xty)  >  0  if  x  ft  ys 

(2.1) 

Hi)  d(x,y)  =  d(ytx ), 
iv )  d(x>y)  £  d(x,z)  +  d(zty). 


The  ordered  pair  (S,d)  is  then  termed  a  metric  space. 


r. 
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If,  for  example,  S  is  the  set  of  real  numbers,  the 
usual  metric  is  d(x,y)  =  | x-y |  .  Properties  (i)  to  (iv) 

are  then  simply  familiar  properties  of  the  absolute  value 
function . 

The  concept  of  a  metric  space  is  usually  too  general 
to  yield  many  practical  results  and  therefore  the  setting 
for  many  problems  in  the  theory  of  approximation  is  a  metric 
space  of  a  particular  kind  called  a  normed  linear  space. 

We  must  first,  however,  define  a  linear  space. 

2.1.2  Linear  Spaces 

Definition  2.2.  Let  S  be  a  set  of  elements 
re,  y ,  a,  .  .  .  for  which  two  types  of  operations  are  possible . 

S  is  termed  a  linear  space  if  any  two  elements  x,y  in 
S  determine  a  unique  element  x+y  in  S  as  their  sum , 
if  each  element  x  in  S  and  each  scalar  a  determines 
a  unique  element  ax  in  S  as  their  scalar  product ,  and 
if  summation  and  scalar  multiplication  satisfy  the  following 
properties : 

i)  x+y  =  y+x , 
ii)  x  +  (y+z)  =  (x+y)  +  z. 

Hi)  There  exists  a  unique  element  0  in  S 
such  that  x+Q  =  x  for  all  x  in  S 
(0  is  called  the  zero  element  of  the  space). 


: 

eqo*xq 

. 

.  ■ 

■ 
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(2.2) 


iv  )  For  each  x  in  S  t  there  exists  a  unique 

inverse  (-x)  in  S  such  that  x  +  (-x)  =  Qs 

v)  a(Qx)  =  (a$)x  for  all  scalars  a,$, 

vi)  a (x+y)  =  ax  +  ay  > 

vii)  (a+$)x  =  ax  +  &Xj 

viii)  lx  =  x. 


We  now  list  some  linear  spaces  which  are  important  to  this 
study . 


Example  2.1.  The  real t  n- dimensional ,  Euclidian 
space j  denoted  by  R n  . 

Rn  consists  of  vectors  x,y,z, . . .  that  are  n-tuples  of 
real  numbers:  x  =  (x^,...,x  )  .  Summation  is  defined  by 
x+y  =  ( x-|_+y ^ i  •  •  •  Jxn+yn )  ,  and  scalar  multiplication  by 
ax  =  (ax^ . . . >axn)  .  The  zero  element  is  0  =  (0,...,0)  , 
and  the  inverse  of  x  is  (-x)  =  (-x^,...,-x  )  . 


Example  2.2.  The  space  of  all  functions 
continuous  in  the  interval  [a_,£>]  >  denoted  by  C\_asb~\  . 

Let  f,g,...  be  members  of  C[a,b]  and  let  x  be  any 
point  in  [a,b]  .  We  define  the  sum  of  f  and  g  by 
(f+g)(x)  =  f(x)  +  g(x)  ,  and  the  scalar  product  of  a  and 
f  by  (af)(x)  =  af(x)  .  The  zero  element  is  that  function 
which  vanishes  identically  in  [a,b]  and  the  inverse  of 
f  is  given  by  (-f)(x)  =  -f(x)  . 


.. 
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Example  2 . 3 .  The  space  of  all  functions  continuous 
on  a  simple  closed  curve  C  ,  denoted  by  W(C)  . 

Example  2.4.  Let  C  be  a  simple  closed  curve  and 
R  the  region  interior  to  C  .  The  space  of  all  complex 
functions  analytic  in  R  and  continuous  on  C  3  denoted  by 
A  (  Rj  C )  . 

Before  leaving  the  concept  of  a  linear  space,  we  shall 
present  a  definition  which  will  be  used  in  subsequent  work. 

Definition  2.3«  bet  S  be  a  linear  space.  A 
subset  M  of  S  is  called  convex  if  xsy  are  any  two 
distinct  elements  of  M  implies  that  qlx  +  (l-a)y  is  also 
a  member  of  M  for  any  real  a  ,  0  <  a  <  1  . 

If  we  picture  M  geometrically,  then  ax  +  (l-a)y 
lies  on  the  line  segment  joining  x  and  y  .  Thus  M  is 
convex  if  it  contains  all  the  points  on  the  line  segment 
joining  any  two  of  its  points.  Circles,  ellipses,  rect¬ 
angles,  and  straight  lines  are  examples  of  convex  sets  in 

R  .  The  figure  "8"  is  not  convex  in  R  . 
n  &  n 

2.1.3  Normed  Linear  Spaces 

Definition  2.4.  A  normed  linear  space  is  a  linear 
space  S  for  which  there  is  defined  a  real-valued  function 
on  elements  x>y3z3...  in  S  called  a  normt  denoted  by 
||x||  ,  and  satisfying  the  following  laws: 
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(2.3) 


i) 

II  ©1  -  0, 

ii) 

Ml  >  o 

if  * 

Hi) 

II  HI  = 

a|  |*| 

for  any  scalar 

iv) 

\\x+y\\  * 

hi  + 

II 2/ II  • 

verified  property  of 

norms  is  given  by 

Theorem 

2.1. 

1*1  - 

II 2/1  *  \x-y\  • 

There  are  usually  an  infinity  of  norms  that  may  be 
introduced  into  linear  spaces  in  order  to  make  them  normed 
linear  spaces.  However,  we  discuss  only  those  relevant 
to  this  study. 


Theorem  2.2. 

defined  by 


The  linear  space  R  with  norm 
e  n 


(2.4) 


||  x\\  =  max  |  x  . 
1  <i<n 


is  a  normed  linear  space. 

Proof:  We  must  verify  that  the  function  given 

in  (2.4)  satisfies  the  properties  listed  in  (2.3).  The 
first  three  properties  are  easily  verified  as  simple 
consequences  of  properties  of  the  absolute  value  function. 
To  prove  (iv)  we  write 


' 


'■a  .  r : : : 

id  ■  •  ©■  .  v 


in 


II x+y II  =  II  (^1+y1»  •  •  •  >xn+yn)ll 

=  max  | x . +y .  | 

l<i<n  1 

=  |xk+yk|  ,  for  some  k 

s  lxkl  +  Ukl 

<  max  | x . |  +  max  | y . 
l<i<n  1  l<i<n  1 

=  Ml  +  llyll  • 

Theorem  2.3.  The  linear  space 

defined  by 

(2.5)  ||  /||  =  max  \  f(x)  | 

a<x<b 


is  a  normed  linear  space. 


Theorem  2.4.  The  linear  space 


defined  by 


(2.6)  || /||  =  max  |  f(x>y) 

(x>y )  zC 


,  l<k<n 


C\_a,b~\  with  norm 


W ( C )  wi th  norm 


is  a  normed  linear  space. 


■ 


15 


Theorem  2.5.  The  linear  space  A(RSC)  with  norm 

defined  by 

(2.7)  ||  /||  =  max  |  f(z)  | 

zeC 


is  a  normed  linear  space. 

Proof:  The  maximum  modulus  need  only  be  taken 

over  C  ,  since  a  complex  function  analytic  in  R  and 
continuous  on  C  assumes  its  maximum  modulus  on  C  ([25], 
pp.  290-291). 

Theorem  2.6.  In  a  normed  linear  space 


(2.8) 


d(f,g)  =  \f-g\ 


defines  a  metric. 


Proof:  We  must  show  that  (2.3)  and  (2.8)  imply 


(2.1). 

i)  d(f,f)  =  |f-f|| 
ii)  d(f,g)  =  )|  f-sl| 
ill)  d(f,g)  =  ||f-g|| 
iv)  d(f,g)  =  ||f-g|| 


=  II 0  •  f|  =  1 0 1  1  f  I  =0  . 

>  0  since  f  i-  g  implies  f-g  t  0  . 

=  1-1. 1  |l  f-g  II  =  lls-fll  =  d(g,f) 

=  ||f-h+h-g|  <  |f-h|  +  |h-g| 


=  d(f,h)  +  d(h,g)  . 


:1oot 

suounlJnoo 


L  it  qml  a  *  i  yon f  0  <  jj  •*  (*  ,  )  i 
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2.1.4  Best  Approximation  From  the  remarks  made  at  the 
beginning  of  this  Section  and  from  Theorem  2.6,  we  see  that 
a  meaningful  way  of  defining  the  closeness  of  two  functions 
f  and  g  in  a  normed  linear  space  is  to  define  it  as  the 
norm  of  their  difference. 


Definition  2.5. 
a  normed  linear  space  S  . 
linear  combinations  of  g 
for •  which 


Let  f,gv  ...,gn  be  e 

A  best  approximation 

...,q  is  an  element 
*  s  n 


lements  of 
to  f  by 


n 


l 

0=1 


c  .q  . 
0*0 


(2.9)  II  f  -  l  OjQA  *  Ilf  "  I  <z  .y  .|| 

0=1  J  J  0  =  1  0  0 


for  every  choice  of  constants  . 

Then  the  fundamental  assertion  of  linear  approximation 
theory  is 

n 

Theorem  2.7.  The  problem  of  finding  min  ||  /  -  £  a  Ag  „.|| 

-  a  .  •  7  J  J 

d  d  1 

has  a  solution . 

A  proof  of  this  important  theorem  can  be  found  in  ([13], 
pp.  137-139). 

Definition  2.6.  For  a  given  f3 g^ , . . . , gn  in  a 
normed  linear  space ,  E  (f)  is  defined  by 

"Pit 

En(f>  =  "in  11/  -  j  Vjl  • 

0  0  J- 


(2.10) 


dx  g|  n 


A  -  ...  .  ,  t  $  'j  ■-  >  \v  :  ./■  '  c  ■  \  j-  j«  -..:  - 
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Theorem  2.8. 


(2.11) 


E1(f)  *  E2(f)  > 


Proof :  (2.11)  is  true  since  linear  combinations 

of  g-j_*«**>gn  are  also  linear  combinations  of  Sq* • • • >§n>Sn+i 

We  have  observed  there  is  always  one  best  approximation 

n 


to  f  by  linear  combinations 
be  more  than  one.  If  M  is 
then  we  have 


l  a.g.  .  However,  there  may 
J-l  J  J 

the  set  of  best  approximations. 


Theorem  2.9.  M  is  a  convex  set. 

A  proof  of  the  theorem  can  be  found  in  ([13],  pp .  140-141). 

In  our  discussion  of  best  approximation,  the  particular 
norm  used  was  not  specified.  However,  since  the  norms 
given  in  Theorems  2. 2-2. 5  are  formulated  in  terms  of  a 
certain  maximum,  the  corresponding  best  approximations  solve 
the  problem  of  minimizing  the  maximum  error,  and  therefore 
we  will  subsequently  refer  to  these  best  approximations  as 
minimax  approximations. 


2 . 2  Some  Properties  of  Harmonic  Functions 

Laplace's  equation  bears  a  close  relationship  to  the 
theory  of  analytic  functions.  In  fact,  we  have 

Theorem  2,10.  The  real  part  of  any  analytic 


function  is  a  harmonic  function;  and  conversely ,  every 


...  s  c\> 

. 

. 

. 
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harmonic  function  is  the  real  part  of  some  analytic  function . 

Proof:  Set  z  =  x  +  iy  and  let  f(z)  =  u(x,y)  + 
iv(x,y)  be  an  analytic  function  of  z  .  Then  u  and  v 
are  related  by  the  Cauchy-Riemann  equations 


(2.12) 


» 


> 


and  are  said  to  be  conjugate  ([25],  p.  83).  Furthermore,  u 

and  v  possess  continuous  partial  derivatives  of  all  orders 

([25],  pp .  185-186).  If  the  first  of  these  equations  is 

differentiated  with  respect  to  x  ,  the  second  with  respect 

to  y,  and  the  two  resulting  equations  added,  then  it  can  be 

seen  that  u(x,y)  =  Re  f(z)  is  a  solution  of  Laplace’s 

equation.  On  the  other  hand,  let  u  be  a  harmonic  function. 

Then  equations  (2.12)  can  be  solved  for  v  ,  and  since  the 

partial  derivatives  u  ,u  ,v  ,v  are  continuous, 
y  x*  y  x’  y  * 

f(z)  =  u(x,y)  +  iv(x,y)  is  an  analytic  function  of  z 
([25],  p.  85). 

A  parallel  argument  yields  a  similar  theorem  for  the 
imaginary  parts  of  analytic  functions. 

Another  very  important  property  of  harmonic  functions 
is  given  in 


Theorem  2.11.  Let  C  be  a  simple  closed  curve 
and  R  the  region  interior  to  C  .  Let  u  be  a  harmonic 


, 

. 

•  i  ■„  .  '  i  ■ 


. 
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function  in  R  and  continuous  on  C  .  Then  there  exist 

two  points  z1  and  z2  on  C  such  that  u(z^)  <  u(z)  <  u(z^) 

for  all  z  in  R  . 

A  proof  of  this  theorem  can  be  found  in  ([25],  pp.  349-350). 

A  crucial  result  of  this  theorem  is 

Theorem  2.12.  Consider  Dirichlet's  problem  for 
Laplace's  equation  (see  Section  1,1).  If  g(x9y)  is  a 
harmonic  function  such  that 

(2.13)  max  | f(xsy)  -  g(x9y) \  <  M 

(x3 y ) eC 

for  some  constant  M  3  then 

(2.14)  max  \u(x9y)  -  g(x9y)  \  <  M  . 

(x9y ) zR 

Proof:  If  u  is  the  solution  of  the  original 

problem,  then  u-g  is  the  solution  of  the  same  problem  with 
boundary  function  f-g  rather  than  f  .  Note  that  f-g 
is  continuous  and  u-g  is  harmonic.  However,  by  Theorem  2.11, 
a  harmonic  function  assumes  its  maximum  value  on  the  boundary 
C  .  That  is , 


max  |u(x,y) 
( x,y ) eR 


g ( x ,y ) |  <  max  |f(x,y)  -  g(x,y)|  <  M  . 
(x,y)eC 


' 
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2 . 3  A  Complete  Set  of  Harmonic  Functions 

So  far  we  have  not  specifically  given  a  set  of  harmonic 
functions  g1,...,gn  with  which  to  approximate  u  in  R 
and  f  on  C  .  We  see,  from  Theorem  2.10,  that  we  need  only 
use  the  real  and/or  imaginary  parts  of  some  analytic  functions. 
Now  zn  ,  n=0,l,...  are  obviously  analytic  functions  and 
hence 

(2.15)  H  =  {1,  Re  zn,  Im  zn,  n=l,2,...} 


is  an  infinite  set  of  harmonic  functions.  They  are  commonly 
called  harmonic  polynomials. 

Once  a  particular  set  of  harmonic  functions  have  been 
chosen,  a  problem  arises  as  to  the  "completeness”  of  this 
set.  Completeness  refers  to  the  possibility  of  arbitrarily 
close  approximation  of  the  boundary  data  by  means  of  the 
boundary  values  of  linear  combinations  of  the  selected  set 
of  harmonic  functions.  A  relevant  result,  due  to  J.L.  Walsh, 
is  that  the  powers,  l,z,z  ,...  are  complete  in  A(R,C) 
([13],  pp .  275-278).  An  important  consequence  of  Walsh’s 
assertion  is 

Theorem  2.13.  The  harmonic  polynomials  are 


complete  in  W(C)  . 


'  '  !. 
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Proof :  Let  u  be  the  solution  of  the  Dirichlet 

problem.  Then  u  is  the  real  part  of  some  analytic 

function  F(z)  ,  by  Theorem  2.10.  By  the  completeness  of 

2 

the  powers  l,z,z  ,  ...  in  A(R,C)  ,  there  exists  a  complex 

polynomial  p  (z)  such  that  max  |F(z)  -  p  (z)|  <  6  for 

zeC 

any  6  >  0  .  We  now  write 


6  >  max  | F ( z )  -  p  ( z  )  | 
zeC  n 

>  | F ( z  )  -  p  ( z ) |  for  any  z  on  C 

>  | Re  (F(z)  -  pn(z) ) I 

>  | f (x,y )  -  Re  pn(z ) | 

Since  the  above  inequality  is  true  for  any  z  on  C  ,  we 
have  that 


6  >  max  | f ( x ,y )  -  Re  p  (z) 
( x ,y  )  eC 


Letting  6  tend  to  zero  proves  the  theorem. 


CHAPTER  III 


MINIMAX  APPROXIMATIONS  ON  A  DISCRETE  POINT  SET 

3 . 1  Introduction 

Given  the  Dirichlet  problem  outlined  in  Chapter  I  (see 
Section  1.1),  we  set  ourselves  the  task  of  selecting  a 
vector  of  parameters  c  =  (c^,...,c  )  to  minimize  the 
expression 

n 

(3.1)  M1  =  max  |f(x,y)  -  l  a.g. (x,y)|  . 

(x,y)eC  j=l  J  3 


Our  discussions  in  Chapter  II  ensured  us  that  such  a  vector 

c  existed.  Furthermore,  if  g.  ,  j=l,...,n  were  chosen 

J 

to  be  the  harmonic  polynomials,  then  for  sufficiently 
n 

large  n  ,  £  c.g.  could  be  accepted  as  a  reasonably 


j=l 


JTJ 


good  approximate  solution  to  the  Dirichlet  problem. 

As  yet,  however,  we  have  not  tackled  the  problem  of 
how  to  determine  the  vector  c  .  As  a  first  step,  we  might 
select  m  points  (x^,y^)  ,  i=l,...,m  (m>n)  from  the 
boundary  C  and  minimize  the  expression 


n 


(3.2) 


M2  = 


max 

l<i<m 


-  i  ajgj (x1>y1> 
J  * 


. 


. 
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accepting  the  minimax  solution  to  (3.2)  as  a  close 
approximation  to  the  minimax  solution  of  (3.1).  This  is 
equivalent  to  solving  the  following  linear  system  of 
equations  in  the  minimax  sense: 


n 


(3.3) 


j=l 


ajgj  (x1,y±) 


=  f (xi,yi )  ,  i=l, .  .  .  ,m 


Obtaining  a  minimax  solution  to  an  overdetermined, 
inconsistent  system  of  linear  equations  is  an  important 
problem  in  its  own  right.  The  discussion  we  now  give 
parallels  that  given  by  Handscomb  ([18],  pp .  73-82).  A 
more  advanced  and  complete  approach  has  been  given  by 
Cheney  ([3],  PP •  28-56). 

3.2  Minimax  Solution  of  an  Overdetermined,  Inconsistent 
System  of  Linear  Equations 

Consider  the  system  of  linear  equations 


n 

( 3 •  *0  I  ai i x i  -  ^i  *  i-l,...,m  . 

j  =  l  J  J 


We  assume  that  the  numbers  a. .  and  b.  are  given  and 

1 J  1 

that  the  unknowns  x.  are  to  be  determined.  (3«*0  is 

J 

frequently  written  in  matrix  notation  as  Ax  =  b  .  The 
system  may  have  exactly  one  solution,  infinitely  many 


'  '  '  .  .  .  : 


. 
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solutions,  or  no  solution,  depending  on  the  data.  If  the 
rank  of  the  coefficient  matrix  A  is  equal  to  the  rank 
of  the  augmented  matrix  ,  then  (3.4)  possesses  a 

solution  and  is  said  to  be  consistent.  If  m>n  and  the 
rank  of  A  does  not  equal  the  rank  of  A^  ,  then  (3.4) 
is  said  to  be  overdetermined  and  inconsistent .  However, 
even  though  (3*4)  may  possess  no  solution,  we  may  still 
try  to  find  a  vector  x  =  (x-^,...,x  )  that  minimizes  the 
expression 

n 

(3.5)  p  =  max  |b.  -  l  a, ,x. |  . 

l<i<m  x  j=l  J 

Such  an  x  is  called  a  minimax  solution  of  (3.4). 

3.2.1  Some  Properties  of  Minimax  Solutions 

Theorem  3.1*  Let  values  a^.^.jb;  be  given  for 
"  1*3  i* 

and  i=l , . , . im>n  .  The  problem  of  finding 


n 

min  max  I  b  .  -  )  a  .  .x  . 

xeR  l<i<m  ^  j  =  l  ^  ^ 

n 


has  a  solution. 

Proof :  For  each  j=l,...,n  ,  let  A^  denote  the 

column  vector  ( alj.  , .  .  .  ,a  .  )T  .  Then  (3.4)  can  be  rewritten 


■ 
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in  the  form 

n 

(3.6)  J  x . A .  =  b  . 

j=l  J  J 

n 

Now  l  y.A.  and  b  are  elements  of  Rm  for  any 
j  =  l  J  J 

y  =  (y1>...>yn)  .  Hence,  by  Theorems  2.2  and  2.7,  there 
exists  a  vector  xeRn  such  that 

(3.7)  II  b  -  l  x  A  ||  <  ||  b  -  l  y  A  || 

j=l  J  J  j=l  J  J 

for  any  yeRn  • 

Theorem  3.2.  The  set  of  minimax  xolutions  to 
(3.4)  is  convex. 

Proof:  Theorem  3.2  is  just  a  rewording  of 

Theorem  2.9. 

Before  proceeding  further,  we  make  the  following 
important  definition. 

Definition  3.1.  A  set  of  vectors  in  R  is  said 
- =: —  j  n 

to  satisfy  Haar’s  condition  if  every  subset  of  n  of  them 
is  linearly  independent . 

We  now  make  the  rather  strong  assumption  that  the  set 
of  row  vectors  making  up  the  matrix  A  satisfy  Haar's 
condition.  Therefore,  any  n  of  the  m  equations  in  (3.*0 


*t  rj  §>.?«.  i  S  k.v 
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possess  a  nonsingular  coefficient  matrix.  Haar’s  condition 
will  be  discussed  in  detail  in  Chapter  IV. 


Theorem  3.3.  Let  x  =  (x^...3x_)  be  a  minimax 

1  J.  YL 


solution  to  (3.4)  and  let  p  =  max  \ b 


n 


n 


1  <i<m 


T  a  .  .x  . 
0=1  ^  0 


Then  I  =  {i  :  I  b.  -  J  a..x.\  =  p}  contains  at  least 

'  I  jlj  13  o' 

(n+1)  elements . 


Proof:  Suppose,  on  the  contrary,  that  I  contains 

r  elements,  where  r  <  n+1  .  Reordering  the  equations  if 
necessary,  we  can  assume  that  the  maximum  absolute  error 
occurs  in  the  first  r  equations.  Thus  we  have 


n 

I 

J-l 


aljXj 


P 


b  -  y  a  ,  x .  =  s  p 
r  r*j  j  rK 

(3.8) 

n 

br+l  “  J,  ar+l , j xj  =  sr+ler+l 


l=1  amj  XJ  smera 
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where 

(3.9)  0  <  ei  <  p  ,  i=r+l,...,m  , 

and 

n 

We  now  show  that  the  rank  of  the  matrix  =  (a^)  > 

j=l,...,n  ;  1=1,. ..,r  ;  is  less  than  r  .  If,  on  the 
other  hand,  the  rank  of  equals  r  ,  then  we  can  find 

a  vector  u  =  (u,,...,u  )  such  that 

(3.1D  aij  uj  “  aiP  .  i“l.  •••.»■  • 

Now  define  a  new  vector  y  =  (y^,...,y  )  by 


(3.10)  s,  =  sgn(b .  -  l  a,  , 

J-l  J 


(3.12) 


y j  —  x j  +  ^  1 , .  .  .  ,  n  , 


where 


P-M-l 

~Mj“  » 


(3.13) 


0  <  X  <  min( 


1)  . 


. 


M1  and  in  (3.13)  are  given  by 


(3.14) 

and 

(3.15) 


Then  for  i=l 


n 


(3.16) 


M-.  =  max  e 
r+l<i<m 


n 

max  I  T  a.  .u .  I  . 
r+l<i<m  j=l  1J  J 


, . . . ,m  we  have 


n 


aijyj 


=  b!  -  I 

j  =  l 


aiJ (XJ 


XUj) 


n 


=  bi  -  I 

j=i 


aijxj 


x  X  a^uj 


/ 

S  ^  P  "  X  S  ^  P  y  l**!^**#^!1  • 


•l»l  - 1 J,  *1,“J 


i=r+l , 


Taking  the  absolute  value  of  (3.16)  we  get 


29 


(3.17) 


n 


I 

j=l 


( 1-X)p  ,  i=l, ...,r  . 


siei 


-  X 


n 

l 

j=i 


a.  .u . 
lj  J 


i=r+l , .  .  .  ,m 


Certainly,  (l-X)p  <  p  and  for  i=r+l,...,m 


siei 


x  X  auuJ! 


5  +  M  I 

j=i 


< 


n 


I  I 

J-l 


aijUj 


(3.18) 


<  ei  +  p  -  M1  ,  by  (3.15) 


*  P  ,  by  (3.1*0. 


This  means  that  y  is  a  better  minimax  approximation  than 
x  ,  which  is  a  contradiction.  Hence,  the  rank  of  A  is 
less  than  r  .  Since  r  <  n  ,  this  contradicts  our  basic 
assumption  that  the  rows  of  A  satisfy  Haar’s  condition. 
Hence,  I  contains  at  least  (n+1)  elements. 
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unique . 


Theorem  3.^.  The  minimax  solution  to  (3.4)  is 


Proof:  Suppose  that  x,y  e  Rn  are  two  distinct 

minimax  solutions  and  let 


n  n 

(3-19)  p  =  max  |b.  -  l  a, .x  |  =  max  |b,  -  £  a, .y,|  . 

l<i<m  1  j=l  1J  J  l<i<m  1  j=l  J 


Now  define 


(3.20) 


v  =  ^(x+y)  . 


Then  v  is  also  a  minimax  solution  by  Theorem  3.2.  For 
those  i  for  which 


(3.21) 


bi  -  X  auvj  ■  sip 


where  s^  is  still  given  by  (3.10),  we  must  have 


n 

I 


n 


(3.22)  (b±  -  ^  aijxj)  +  (bi  “  aijyj)  ’  2sip 


Therefore,  we  obtain 
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(3.23)  2p  <  |b±  -  ^  aijXj|  +  |b1  -  ^  aijyj|  . 


Neither  term  on  the  right  of  the  above  inequality  can  be 
greater  than  p  .  Therefore, 


n 


n 


(3.24)  2p  =  |b.  -  a^.x.l  +  |b±  -  ^  a±jyj 


j=l 


j  =  l 


By  the  same  reasoning,  we  obtain 


n 


n 


(3.25)  P  *  |bj_  -  l  a  x  |  =  |bi  -  I  a±.y 


j=l 


j=l 


* 


Hence , 


(3.26) 


n 

-  I 


n 


1  -  ji-L  aUXj  '  bi  '  J],  aijyj  * 


for  if  the  two  terms  were  of  opposite  sign,  this  would 
contradict  (3*22).  Therefore, 


n 


E  au(xj  -  yf  =  0  • 


j=i 


(3.27) 


- 


■  ■  ■:.!  .  ;  • 
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By  Theorem  3.3,  there  are  (n+1)  such  i  for  which  (3.27) 
holds.  Choose  any  n  of  them.  The  above  equation  is 
then  a  system  of  n  equations  in  n  unknowns.  Our  basic 
assumption  ensures  that  the  coefficient  matrix  of  this 
system  is  nonsingular.  Hence,  (3.27)  has  only  the  zero 
solution.  That  is. 


(3.28) 


xj  -  yj  *  0  *  . 


and  therefore  x=y  as  was  to  be  shown. 

Now  that  we  have  some  insight  into  the  existence, 
characterization,  and  uniqueness  of  the  minimax  solution 
to  ( 3 .  *0  ,  we  must  approach  the  problem  of  actually 
computing  it.  We  first  turn  our  attention  to  the  special 
case  where  the  number  of  equations  is  one  greater  than 
the  number  of  unknowns. 

3.2.2  The  Minimax  Solution  of  (n+1)  Equations  in 
n  Unknowns  From  Theorem  3*3,  we  know  the  minimax  solution 
yields  the  same  absolute  error  in  each  equation.  Call 
this  absolute  error  p  .  Using  the  notation  s^  *  ±1  for 

Vi 

the  sign  of  the  error  in  the  i  equation,  we  have  that 
the  minimax  solution  x  satisfies 


n 


(3.29) 


-  I 


aijXJ 


sip 


i»l, . . . ,n+l 
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For  a  particular  choice  of  signs,  the  above  system  can  be 
rewritten  as 


(3.30) 


i  =  l , . . , ,n+l  , 


which  is  simply  a  system  of  (n+1)  equations  in  the  (n+1) 
unknowns  x^,.,.,xn,p  .  Let  us  define 


(3.3D  T,  -  (-1) 


a 


11 


a 


n+1 , 1 


a 


In 


ai-l,l  al-l,n 


al+l,l  "•  al+l,n 


a 


n+1  ,n 


,  1 - 1 i i • i | n+ 1 


Note  that  by  Haar's  condition,  t  0  for  all  i  .  We 
then  have 

Theorem  3i5.  The  eigne  of  are  *^e 

eame  ae  or  opposite  to  the  oorreeponding  eigne  of 

T1  • 
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Proof:  The  coefficient  matrix  of  (3*30)  is 

given  by 


'll 


'In 


n+1,1  * 


a  .  -i  s  , , 
n+l,n  n+1 


Solving  for  p  by  Cramer’s  rule,  expanding  determinants 
down  the  last  column,  and  utilizing  (3. 3D,  we  get 


(3.32) 


s  .  T . 

l  l 


We  are  trying  to  determine  s^  so  that  p  is  minimized. 
Hence,  we  must  try  to  make  the  denonimator  in  (3.32)  as 
large  in  absolute  value  as  possible.  This  will  happen  if 


(3.33) 


s^  =  sgn  T^  ,  i=l,.,.,n+l  , 


or 


(3.34) 


s^  =  -  sgn  T^  ,  i=l,.,.,n+l  , 


which  was  to  be  shown. 


■ 
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Now  that  the  signs  s^  have  been  independently 
determined,  it  remains  to  be  shown  that  the  coefficient 
matrix  of  (3.30)  is  nonsingular. 

Theorem  3.6 


(3.35) 


D  = 


a 


11 


a 


1  n 


s 


a 


n+1 ,  1 


a  .  i  s  .  i 

n+l^n  n+1 


ft  0 


Proof:  We  have  already  remarked  that  each 

T.  ^  0  .  To  prove  the  assertion,  we  expand  the  determin¬ 
ant  down  the  last  column,  obtaining 


n+1 

(3-36)  D  =  l  s.T. 

1-1 


But  by  Theorem  3-5,  we  have  chosen  s^  according  to  (3.33) 
or  (3.3*0,  and  in  one  case,  we  have  a  sum  of  all  positive 
terms,  and  in  the  other,  a  sum  of  all  negative  terms.  In 
either  case, 

(3.37)  D  jt  0  , 


and  the  theorem  is  established. 
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To  summarize,  we  can  obtain  the  minimax  solution  of 
(n+1)  equations  in  n  unknowns  by  first  determining  the 
numbers  s.  according  to  (3*33)  or  (3.3*0  and  then  solving 
the  linear  system  (3-30)  for  x19...,x  ,p  . 

We  now  return  to  the  original  problem  of  determining 
the  minimax  solution  to  (3**0.  The  results  just  presented 
are  important  because  the  minimax  solution  of  ( 3 .  *0  is 
also  the  minimax  solution  to  some  subsystem  of  ( 3 .  *0 
consisting  of  just  (n+1)  equations  (see  Theorem  3*3). 

The  main  problem  now,  then,  is  to  determine  which  subset 
of  (n+1)  equations  is  pertinent.  This  subset  can  be 
determined  by  the  exchange  algorithm. 

3.2.3  The  Exchange  Algorithm  The  basic  idea  of  this 
algorithm  is  to  calculate  the  minimax  solutions  of  a 
succession  of  subsystems,  each  consisting  of  (n+1) 
equations  from  ( 3 .  *0  •  By  Theorem  3-3,  the  minimax  solution 
of  one  of  these  subsystems  is  the  required  minimax  solution. 
Since  there  are  only  a  finite  number  of  possible  sub¬ 
systems,  the  exchange  algorithm  will  converge  in  a  finite 
number  of  steps,  providing  we  can  show  that  it  never 
reconsiders  a  previous  subsystem. 

3.2.3* 1  Theory  of  the  Exchange  Algorithm 
Initially,  we  choose  any  (n+1)  equations  from  ( 3 .  *0  • 
Without  loss  of  generality,  we  may  assume  that  we  have 
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selected  the  first  (n+1)  .  Then  form  the  and  s^ 

according  to 


(3-38) 


T.  =  (-1) 


'll 


i-1,1 


i+1 , 1 


n+1, 1 


In 


i-l,n 


i+1  ,n 


n+1  ,n 


,  i=l, . . . ,n+l  , 


and 


(3.39) 


s^  =  sgn  ,  i=l,...,n+l  . 


The  values  of  the  (n+1)  unknowns  xi>**’>xn>P  are  obtained, 
by  solving 


n 


(3.^0)  l  a. .x  +  s.p  =  b.  ,  i=l, . . . ,n+l  . 

j=i  ij  j  1  1 


We  may  at  this  point  assume  that  p  >  0  .  If  (3.^0)  gives 
p  <  0  ,  we  can  simply  replace  each  s^  by  -s ^  and  each 
T^  by  -T^  .  We  now  compute  the  discrepancies  in  each 
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equation  from 


n 


(3. 'll) 


f  •  .  •  ,  m  j 


and  determine  max  p.  .  A  test  is  made  to  see  if 

l<i<m  1 

max  p.  >  p  .  If  not,  then  x  is  our  required  minimax 
l<i <m  1 

solution  and  p  is  the  associated  minimax  error.  On  the 
other  hand,  if  the  above  condition  is  true,  then  we  may 
assume  that  max  p^  =  Pn+2  >  so  that 


I<i<m 


n 


where,  as  usual,  sn+2  is  the  s^Sn  the  error  in  the 
(n+2)nd  equation. 


Theorem  3»7.  If  (3.42)  holds ,  then  the  minimax 


error  p'  of  some  subsystem  consisting  of  (n+1)  equations 


selected  from  equations  l3...3n+2  is  greater  than  p  . 

The  method  of  proof  will  be  to  explicitly  give  the 
indices  of  the  equations  comprising  the  relevant  subsystem. 


Proof:  We  define  a  set  of  row  vectors  by 


(3. 43) 


Ai  (ail’ - ‘ ,ain)  ’  1’1 


,  .  .  .  , 


n+2  . 


•wxis  x£mi 


. 


.. 


1 
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Now  consider  the  determinants  given  by 


(3.44)  D  = 

J 


L11 


'In 


an+l ,  j  an+l,l  •“  an+l,n 


j  —  1 ,  • • . ,n 


Obviously,  Dj  =  0  for  all  j  since  the  1st  and  (j+l)st 
columns  of  the  determinant  are  identical.  If  we  expand 
D.  down  the  first  column,  we  obtain 


(3.45) 


n+1 

°  =  ±l1  aijTi  *  j=1, 


> 


which  may  be  rewritten  as 


n+1 

(3.46)  I  T.A,  =  0  . 

i=l  1  1 


Furthermore,  since 


1  ,  we  may  rewrite  (3.46)  as 


(3.47) 


n+1 

I  TiS.CsjAO  =  0  . 
i-1  1  11 


If  we  consider  (3.47)  as  a  linear  combination  of  the  vectors 
s^A^  ,  then,  by  (3.39),  all  the  coefficients  are  positive. 
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Since  no  T.  =  0  ,  we  may  solve  for  any  A  from  (3.46), 

1  K 

getting 


(3.48) 


n+1  T. 

=  I  "  T~  ^i  * 
K  i=l  1 

i^k 


Now  by  Haar's  condition,  we  know  that  the  vectors  A^,...^ 
span  Rn  .  Therefore,  there  exist  numbers  Ai  such  that 


(3.49) 


A 


n+2 


n+1 


i  =  l 


xiAi 


Hence , 


n+1 


0  An+2  "  l 


1  =  1 


XiAi 


=  A 


n+2 


XkAk 


n+1 

l  \\ 

i  =  l  K  K 
i^k 


(3.50) 


=  A 


n+2 


n+1  T 

Ak  ^  ”  T~  Ai 

K  i=l  xk  1 

i^k 


n+1 

I  V* 

i  =  l  1  ' 
i^k 


n+1  T, 

An+2  +  (Ak  \  "  Xi)Ai 
i^k 


n+1  T 

+  £  sn+2si^XkT-* 

i=l  n  d  1  K  k 

i?*k 


sn+2An+2 


Xi  )siAi  * 


. 


■ 


n 


The  last  line  comes  from  multiplying  the  previous  one  by 


2 

sn+2  and  aSain  noting  that  s^  =  1  .  The  above  equation 


expresses  the  fact  that  we  have  found  numbers  such  that  a 
linear  combination  of  s^A^  ,  i=l,...,n+2;  i^k  with  these 
numbers  as  coefficients  is  equal  to  zero.  Referring  to 
the  discussion  after  (3.47),  we  must  have  all  of  these 
coefficients  positive.  That  is,  we  require 


(3.51) 


, .  .  .  , 


n+1;  i^k 


Therefore,  we  must  have 


(3.52) 


n+1; 


Prom  (3-39),  we  have  >  0  .  Hence 


(3.53) 


>  sn+2  Ti 


i  =  l 


, .  .  .  , 


n+1;  i^k 


and  therefore  k  is  the  index  satisfying 


(3.54) 
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The  index  k  thus  chosen  is  unique,  for  otherwise  one  of 
the  coefficients  in  (3.50)  would  vanish,  contradicting 
Haar’s  condition.  It  is  now  asserted  that  the  minimax 
error  p'  associated  with  equations  1, . . . ,k-l,k+l, . . . ,n+2 
is  greater  than  the  minimax  error  p  associated  with 
equations  l,..,,n+l  .  Let  x'  be  the  minimax  solution 
of  equations  1, . . , ,k-l,k+l, . . . ,n+2  .  Then 


n  n 


i  =  l 


» •  •  •  i 


n+1  , 


and 


n  n 


i  =  l 


i  •  •  '  > 


n+1;  i^k  . 


Hence,  x  ^  x'  ,  or 


(3.57) 


x  -  x’  ?  0  . 


■  *  K: 

r 

. 


. 
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Now  we  can  write 


n 


n 


si  .1,  ai.i(x.i-x.i)  =  si(bi  "  X  “uV 


j  =  l 


J-l 


(3.58) 


n 


-  s 


i(bi  -  aijxJ 


=  p  -  p '  ,  i=l,...,n+l;  i^k  . 


Now  by  Haar's  condition  and  (3*57),  we  have 


(3.59) 


p  -  p*  ¥■  0  . 


We  also  have  that 


n+2 


n 

l 

j  =  l 


(xt-x4)  = 


n+2 ,  j  j  j 


sn+2(bn+2  "  £ 


n 

I 

j-i 


(3.60) 


"  sn+2(bn+2  “  ^  an+2,jxj) 


>  P  -  P*  • 


. 


. 
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We  are  trying  to  show  that  p  -  p*  <  0  .  We  have  already 
shown  that  p  -  p*  ?  0  .  Suppose,  then,  that  p  -  p’  >0  . 
Then  by  (3.58)  and  (3.60),  we  can  write 

n 

(xj-Xj)  >  0  ,  i=l,...,n+2;  i^k  . 

From  (3-50)  and  (3.54),  we  know  there  exist  numbers  ui 
such  that 


(3.6l)  s,  I  a,, 

j=l  J 


(3.62) 


n+2 

I 

i=l 

i^k 


yisiAi 


=  0 


or  equivalently. 


(3.63) 


n+2 

I 

i  =  l 
i^k 


yisiai j 


=  0 


j  = 1 , . * .  ,n 


Since  the  >  0  ,  we  can  rewrite  (3.61)  as 


Vi  aij(xj-xJ)  >  0  • 


(3-64) 


> 


,n+2;  i^k  . 


'.7  Si.  tV- 
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For  any  numbers  y ^  ,  (3.63)  can  be  altered  to  read 


(3.65) 


n+2 

yj  Jx  nsiAi 

i^k 


—  0  ,  J  1 » •  •  •  j  n 


Summing  the  last  set  of  equations  over  j  ,  we  obtain 


(3.66) 


n+2 

l  ^isi 

i=l  1  1 

i^k 


n 

l 


yjau 


0  . 


The  above  equation  is  true  for  any  set  of  numbers 
In  particular,  then,  it  must  be  true  for  y^  -  xj 
j=l,...,n  ,  and  thus  we  get 


t 


(3.67) 


n+2  n 


l-±  Vi  aij(xj-xj 


)  =  0  . 


i 
i^k 


If  we  sum  (3.64)  over  all  permissable  indices  of  i  ,  we  get 


n+2  n 

J,  Vi  aij(xj-xj)  * 0  • 

i/k 


(3.68) 
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which  contradicts  (3-67).  Therefore,  p  -  p  ’  <0  and  thus 
p’  >  p  as  was  to  be  shown. 

We  can  now  state  that  the  exchange  algorithm  converges 
in  a  finite  number  of  steps.  From  the  previous  theorem, 
we  know  that  the  minimax  error  strictly  increases  during 
each  cycle  of  the  algorithm.  Hence,  the  algorithm  never 
returns  to  a  previously  considered  subsystem.  Since  there 
are  only  a  finite  number  of  such  subsystems,  the  algorithm 
must  eventually  reach  the  critical  one  whose  existence  is 
guaranteed  by  Theorem  3.3. 

3. 2. 3. 2  The  Exchange  Algorithm  on  a  Computer 
Our  approach  in  this  section  will  be  to  present  the 
exchange  algorithm  step  by  step  with  explanations  inserted 
where  necessary. 

From  the  last  section,  we  know  that  we  will  always 
be  considering  some  subset  of  (n+1)  equations.  This 
can  be  most  easily  handled  by  an  (n+1)  -  component  vector 
I  =  (i1 ,i2 , . . . ,in+1 }  which  gives  the  indices  of  the 
(n+1)  equations  presently  under  consideration.  We 
suppose  that  the  matrix  A  ,  the  vector  b  ,  and  a  vector 
I  have  been  given.  The  exchange  algorithm  then  proceeds 


as  follows: 


gnliub  sssBQTonl  v^^oi'de  iotiq  xi:, minim  y  t  •£. .."  v  cm 
nevsn  mri^l'iosla  sritf  « sonsH  *cU  ••.  •  ■  >v  fk  bs 

9*1  ©dd  sonlS  .ms^e^cidws  fceisblanoo  ^XauoX'.  .  .:•>  -  . 

raridX'iogXa  ©rid  t  amsda^edus  rlous  lo  ‘xedourn  edlnll  b  ^Xno  srtfi 
2 t  aonsdalxs  ssorfw  9fto  iBoictXio 

•  ... 

bsdiesnl  anoXd£n£..qx9  rid.iv/  qsds  v-‘  ^  iviltlnc;..£s  ixt* 

‘ioJosv  dnanoqrrjoo  -  (1+n) 

.noIdB'isbXenoo  'u  ..jj  \  l  t  n  .L, 

; 
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1) 


Tk  '  (-d) 


3  •  -i  •  •  •  cl  • 

11*1  1l>n 


i  1 


1k+l»1 


Vl>n 


ik+l,n 


ai  i  *  •  •  ai  n 

■Ln+1» 


n+1 


,  k=l, . . . ,n+l  . 


Calculating  the  determinant  is  a  standard  problem  in  linear 
algebra.  For  this  and  all  other  problems  in  linear  algebra 
we  have  adapted  a  set  of  programs  given  by  Forsythe  and 
Moler  ([15],  pp .  68-72).  Their  main  computational  tool 
is  Gaussian  elimination,  with  partial  pivoting  and  iterative 
improvement.  For  completeness,  these  programs  have  been 
listed  in  the  Appendix. 


2)  sk  <-  sgn  Tk  ,  k=l ,  .  .  .  ,n+l  . 


—  — 

dll  dl,n+l 

—  — 

3  •  -i  •  •  *  3  •  S 

i1#l  i1#n  1 

•  • 

•  • 

•  • 

•  •  • 

•  •  • 

•  •  • 

dn+l , 1  dn+l ,n+l 

ai  i  • • •  ai  n  sn+l 

n+1 * x  1n+l» 
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r 

i — i 

L^_ 

dll  •••  dl,n+l 

• 

• 

•  0 

•  * 

• 

# 

0 

•  0 

0 

X 

n 

dn+l,l  dn+l,n+l 

b. 

1n+l 

P 


The  accumulation  of  inner  products  is  carried  out  in  double 
precision . 


n 

5)  Pi  *  b±  -  I  aij xj  .  . 

6)  | p J  =  max  | p, |  . 

l<i <m  1 

7)  We  test:  Is  Ip  I  >  p  ?  If  the  answer  is  no.  then 
[x^,...,x  ,p]  gives  the  minimax  solution  and  the  associated 
minimax  error.  If  the  answer  is  yes,  we  have  found  it 
necessary  to  make  a  further  test  to  see  if  a  e  I  .  If 

so,  then  we  exit  with  the  vector  [x^,...,x  ,p]  as  before. 
If  not,  then  we  proceed  to  the  next  step. 

8)  u  sgn  p^  . 
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A  word  of  explanation  is  necessary  to  explain  the  following 
steps  in  the  algorithm.  From  the  proof  of  Theorem  3.7, 
we  see  that  we  must  express  s  A  as  a  linear  combination 

of  s]_Aj,  •  •  .  >s  ^A^  To  do  this,  set 

1  n+1 


(3.69) 


n+1 


a 


1 


k=l 


XkAi, 


The  above  equation  can  be  written  in  matrix  form 


!>]_»•  •  •  >An+i] 


(3.70) 


cl#  -i  •••  cl.  S-i 

h-1  ii>n  1 


'nil*1  ’ 


a .  s  , , 

1n+l,n  n+1 


—  C  a  ,  •  •  •  , a  „ ,  u  J  . 
oil*  5  an* 


Therefore,  the  next  step  in  the  algorithm  is  given  by 


[a-l,  . . .  ,xn+1]  -  [ aal >  •  •  •  » aan  1  u ^ 


dll  **•  dl,n+l 


n+1 , 1  •**  n+1, n+1 


- 
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We  must  now  determine  k  from  (3*5*0.  From  the  definition 
of  an  inverse,  we  have 


dll  **•  dl,n+l 

•  • 

3.*  -|  •  •  •  S-. 

i1,l  i14n  1 

•  •  • 

•  • 

dn+l,l  ***  dn+l ,n+l 

•  •  • 

•  •  • 

3. .  i  ...  a.  s.t 

dn+l * 1  in+l,n  n+1 

(3.7D 


1  0  ...  0 

0  1  ...  0 

•  •  • 

•  •  • 

•  e  • 


0 


0 


1 


and  hence. 


(3.72) 


n+1 

I  ^n+l  k^i  =  0  • 
k=!  n  ljK  1k 


Comparing  (3*72)  with  (3.^6),  we  see  that  after  the  initial 
inverse  has  been  calculated,  the  numbers  ^n  +  ]_  ^  may 
serve  as  the  numbers  T  .  From  (3.69)  we  have 


• 

and  therefore,  the  next  step  in  the  algorithm  can  be  written 


Xk  Ai 

10)  u  -3 -  =  max  u  -= -  . 

n+l,k  l<i<n+l  an+l,i 

The  only  thing  left  to  do  is  to  update  the  inverse.  We 
are  changing  the  k  row  of  a  coefficient  matrix  by 
replacing  it  by  a  known  linear  combination  of  rows,  the 
coefficients  being  the  X^  of  (3*69).  The  effect  on  the 
inverse  is  given  by  Handscomb  ([18],  pp .  79-80). 

dij/Xj  »  J=k;  i  =  1»* *  * ,n+1  ‘ 

dij  dik  *  i  =  1>-  •  •  >n+1  ' 


12) 


a  . 


The  algorithm  now  returns  to  step  4. 

A  program  to  compute  the  minimax  solution  of  (3.4) 
is  given  in  the  Appendix. 
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CHAPTER  IV 


HAAR’S  CONDITION 


4 . 1  An  Example 

In  the  previous  chapter  we  presented  a  detailed 
exposition  of  the  exchange  algorithm  for  computing  the 
minimax  solution  of  an  overdetermined,  inconsistent  system 
of  linear  equations  Ax  =  b  (see  Section  3.4).  The 
theoretical  basis  of  this  algorithm  relied  heavily  on  the 
assumption  that  the  row  vectors  comprising  the  coefficient 
matrix  A  satisfied  Haar's  condition  (see  Definition  3.1) • 
To  see  the  difficulties  arising  when  Haar's  condition  is 
violated,  consider  the  problem  of  determining  the  minimax 
solution  of  the  following  system  of  linear  equations 


2xx  +  x2  =  6.9 


3x±  +  x2 


7.2 


3.0 


(4.1) 


x 


1 


X 


2 


1.0 


2x1  +  4x2  =  11.1 


x-^  +  2x2 


7.0 


53 


We  note  that  Haar's  condition  is  violated  since  the  vectors 
(2,4)  and  (1,2)  corresponding  to  the  fifth  and  sixth 
equations  are  linearly  dependent.  At  each  stage,  the 
exchange  algorithm  considers  subsystems  of  three  equations. 
As  in  Section  3. 2. 3. 2,  we  will  let  I  be  a  3-component 
vector  giving  the  indices  of  the  equations  presently  under 
consideration.  Initially,  suppose  I  =  [1,2,3]  .  We 
present  the  computations  of  the  exchange  algorithm  given 
in  Section  3 . 2 . 3 . 2 . 


Step  2 


1)  Tx 


(-1) 


3  1 
1  1 


2)  s1  =  sgn  T±  =  -1,  s 2 


3) 


D 


. 

' 
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4) 

X1 

0 

1 

2 

1 

"2 

6.9 

2.1 

x2 

= 

1 

2 

3 

5 

"? 

7.2 

= 

1.8 

P 

1 

"2 

1 

"4" 

1 

3.0 

-0.9 

5)  p1  =  0.9  ,  p2  =  -0.9  ,  p^  =  -0.9  , 


Pi,  =  0.7  ,  P5  =  -0.3  ,  P6  =  1.3  . 


6 )  | p  r  |  =  max  | p .  |  ;  a  =  6  . 

°  l<i<6  1 


7)  | p6 |  >  P  .  6*1  =  [1,2,3]  . 


8)  u  =  sgn  p6  =  1  . 


9) 

[X1,X2,X3]  =  [1,2,1] 


1  3  5 

2  “TT  ¥ 

111 
"2  T\  ? 


rl  _3  9n 
L  2  >  17*  Tf-* 


10) 


3 


k  =  3  . 


u 

“d 


max  u  ^ — 
l<i<3  a3i 


12)  I  =  [1,2,6]  . 

Step  2 


4) 

xi 

1 

9 

1 

3 

2 

9 

6.9 

29 

IF 

x2 

= 

2 

9 

1 

3 

5 

9 

7.2 

= 

136 

45 

P 

5 

"9 

1 

3 

1 

9 

7.0 

59 

_”9°_ 

5) 


'1 


59.  .  =  _59 

90  J  p2  90 


147 
90  * 


P  4 


217 
90  , 


^  =  379.  n  =  59 
p5  90  5  p6  "90  * 


6 )  |  p  j- 1  =  max  |  p  . 

5  l<i<6  1 


a  =  5  . 
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7)  |  p5  |  >  p  .  5  t  I  . 


8)  u  =  sgnp<-=-l. 


9) 


[X1SA2,A3]  =  [2,4, -1] 


112 

9  3  "9 

2  _1  5 

9  3  9 

5  11 

9  3  9 


10) 


u  A1 


max 

l<i<3 


k  =  1  or  k  =  2  . 


We  have  shown  in  Chapter  III  the  maximum  of  the  above  ratios 
would  be  unique  if  Haar's  condition  were  satisfied.  We  are 
presented  with  two  choices  above,  and  we  will  consider 


both  of  them. 


Case  1  (k  =  2) 


11) 


D  = 


2  11 
3  "3  3 


1  2 
3  3 


12)  I  =  [1,5,6]  . 

Step  3 


4) 

-  — 

X1 

2 

3 

1 

3 

1 

3 

6.9 

97 

30 

x2 

= 

1 

”3 

1 

3 

0 

11.1 

= 

14 

10 

. 

P 

_  _ 

0 

1 

“3 

2 

3_ 

7.0 

29 

30 

_  _1 

5) 


29 

_  39 

49 

30  9 

p2  10  9 

p3 

30 

25 

p6  = 

29 

30  9 

p  = 

p5  30  5 

30 

6 )  I  p0  I  =  max  I p .  I  ;  a  =  2  . 
^  l<i<6  1 


58 


7)  |  p2 I  >  P  •  2  t  I  . 


8)  u  =  sgn  p2  s  -1  . 


9) 


2  1 

3  "3 


[A1,A2,A3]  =  [3, 1,-1] 


1  1 
3  3 


0 


1 

3 


=  [1  i] 

l3>  3J  3J 


10) 


u  A1  =  5/3 


u  A  2 

d32 


-3 


-  1/2  . 


Our  computations  break  down  at  this  point  because  of  the 
division  by  zero.  Let  us  therefore  consider  the  second 
alternative  at  stage  10  in  step  2. 


59 


Case  2  ( k  =  1 ) 


ID 


D  = 


_1  2 
15  5 

__2  1 

15  5 


1 


1 

3 


2 

3 


12)  I  =  [5,2,6]  . 


We  now  carry  out  the  computations  in  Step  3.  However,  the 
computations  will  again  break  down  at  stage  10,  for  in 
calculating  u  X^/d^2  >  we  will  once  more  encounter  a  zero 
divisor.  Thus  the  exchange  algorithm  has  failed  in  this 
example  with  initial  choice  of  I  given  by  I  =  [1,2,3]  . 

There  are  two  possible  ways  of  working  out  of  this 
dilemma.  First,  we  can  try  to  show  in  the  context  of 
the  Dirichlet  problem  that  Haar's  condition  is  automatically 
satisfied.  If  this  approach  fails,  we  can  try  to  modify 
the  exchange  algorithm  so  that  it  will  work  even  though 
Haar’s  condition  is  violated.  Let  us  examine  the  first  of 


these  alternatives. 


60 


4 . 2  Unisolvent  Functions 

Consider  a  set  of  functions  defined  on  a 

point  set  X  .  We  then  have 

Definition  4.1.  The  set  of  functions 

J.  VI 

is  said  to  be  unisolvent  if  the  determinant 


(4.2) 


G  = 


g 1 ( x i)  ...  g  ^ ( x  j ) 


g  n(x  )  ...  g  (x  ) 

*-7  7  ft  &  v?  y) 


n  n 


is  non-zero  for  any  n  distinct  points  x^a...tx  in  X  . 

Suppose  now  that  we  have  a  set  of  unisolvent  functions 
gl,***,gn  on  a  Point  set  x  •  Then  the  coefficient  matrix 
of  the  linear  system 


n 

(4.3)  l  c  g  (x. )  =  f (x. )  ,  i=l,...,m>n 

j  =  l  J  J 


obviously  satisfies  Haar’s  condition. 

In  the  context  of  Dirichlet’s  problem,  the  linear 


system  is 


■ 


■ 


6l 


(4.4) 


n 

I 

j=l 


(x. ,y . )  =  f (x. 


c  .  g  .  x .  ,  y 
J6J  i  4  *  i 


>y± ) 


1=1, . . . ,m>n  , 


where  the  m  points  are  selected  from  the  contour  C  . 
As  a  consequence  of  the  following  assertion,  we  should 
always  approximate  f  by  at  least  three  functions. 

Theorem  4.1.  No  set  of  two  functions  g^tg^ 
are  unisolvent  on  a  simple  closed  curve  C  . 

Proof;  Suppose  the  determinant 


(4.5)  G  = 


gl(xl,yl) 

g2(xi>yi) 

gl(x2,y2) 

g2(x2»y2) 

is  non-zero  for  a  certain  pair  of  points  (x^y^)  and 
(x2,y2)  .  By  a  continuous  motion  around  the  curve  C  , 
these  two  points  can  be  interchanged  without  becoming 
coincident.  In  effect,  we  have  then  interchanged  two 
rows  in  the  determinant  in  (4.5).  Hence  G  has  changed 
sign.  Therefore,  at  some  intermediate  position  the 
determinant  must  vanish,  which  was  to  be  shown. 

Let  us  now  turn  to  the  case  where  g1,*.*Jgn  are  the 
harmonic  polynomials.  We  must  then  examine  the  behavior 
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of  the  following  determinant: 


(4.6) 


1  Re  z 


1 


Im  z 


1 


Re  z^ 


Im  z^ 


G 


n 


1  Re  z 


2n+l 


Im  z 


2n+l 


where  zi>''’>z2n+l  are  distinct  points  from  C  .  To 
show  a  simple  example  of  the  vanishing  of  the  above 
determinant,  suppose  n  =  1  .  Then  (4.6)  reduces  to 


1  xi 


(4.7) 


G1  “  1  x2  y  2 


where  zk  =  xk  +  *  k=l,2,3  .  G^  =  0  if  the  points 

(xl>y^)j  (x2,y2),  (x^,y^)  lie  on  any  straight  line.  There¬ 
fore,  if  C  contains  a  straight  line  segment,  the  functions 
1,  Re  z,  Im  z  are  not  unisolvent  on  C  . 

The  condition  for  the  vanishing  of  Gn  is  that  there 
exist  some  linear  combination  of  the  functions 
1,  Re  z,  Im  z,  ...,  Re  zn,  Im  zn  ,  with  coefficients  not 
all  zero,  which  vanishes  at  the  points  zi**‘->z2n+l  * 
Therefore,  the  condition  for  the  nonvanishing  of  Gn  is 


. 
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that  the  points  zi >  •  ‘  > z2n+l  not  lae  on  any  curve  of 

the  form  H  (z)  =  0  >  where  H  (z)  is  a  linear  combination 
of  the  first  2n+l  harmonic  polynomials.  If  we  wish  this 
to  be  true  for  any  2n+l  points  on  C  ,  then  we  must  have 
that  C  does  not  intersect  any  curve  of  the  form  H  (z)  =  6 
at  more  than  2n  points . 

For  an  arbitrary  simple  closed  curve  C  ,  this  is  a 
difficult  condition  to  check.  Now  H  (z)  =  0  is  an 
algebraic  curve  of  degree  n  .  An  algebraic  curve  of  degree 
n  can  intersect  an  algebraic  curve  of  degree  2  in  at 
most  2n  distinct  points.  Hence,  if  C  is  an  ellipse  or 
a  circle,  the  harmonic  polynomials  are  unisolvent  on  C  . 

A  more  extensive  discussion  of  these  matters  can  be 
found  in  Curtiss  [9,10]. 

4 . 3  Cheney’s  Perturbation  Method 

The  results  of  the  previous  section  are  rather  dis¬ 
couraging  in  that  the  only  interesting  boundaries  on  which 
we  can  be  sure  that  the  harmonic  polynomials  are  unisolvent 
are  circles  and  ellipses.  However,  an  alternative  way  of 
handling  the  problem  is  a  perturbation  technique  suggested 
by  Cheney  ([13],  p.  51).  The  method  is  suggested  by  the 
following  assertion. 

Theorem  4.2.  If  for  a  square  matrix  A  =  ( a j.  .)  , 
-  I'O 

\a I  =  0  ,  then  for  sufficiently  small  e  ^  0  ,  |A-el|  ^  0  . 
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Proof:  Suppose  that  |A  — A I |  =  0  .  Then  A  is  an  eigen¬ 

value  of  the  matrix  (a^)  >  i,j=l,...,n  •  There  are  at 
most  n  such  eigenvalues,  say  A1,...,An  .  If  e  is 

chosen  such  that  0  <  |e|  <  max  |A.|  ,  A.  ^0  ,  then  the 

l<j  <n  J  J 

theorem  is  established. 

To  show  how  this  method  works,  we  again  consider  the 
problem  of  finding  a  minimax  solution  of  (4.1).  Since  the 
vectors  comprising  the  fifth  and  sixth  equations  are 
dependent.  Theorem  4.1  suggests  that  we  first  determine 
the  minimax  solution  to 


2xl 

+ 

X2 

6.9 

1 — 1 

X 

on 

+ 

x2 

7.2 

xi 

+ 

x^  = 

o 

• 

on 

xi 

- 

x2 

1.0 

( 2+e 

+ 

4x^  = 

i — 1 

• 

i — 1 
i — 1 

X1 

+ 

(2+e)x2  = 

7.0 

where  0  <  |e|  <  4  .  We  will  follow  the  same  format  as 
in  Section  4.1.  Initially,  we  again  choose  I  =  [1,2,3]. 


■ 


6  5 


Step  1 

1)  T 1  =  -2,  T2  =  1,  T3  =  1  . 


2)  si  =  “1>  =  1,  =  1  . 


3) 

2 

1 

-1 

-1 

0 

1 

2 

1 

2 

D  = 

3 

1 

1 

= 

1 

2 

3 

T 

5 

¥ 

1 

1 

1 

1 

2 

1 

1 

? 

~xl~ 

"6.9" 

2.1 

X2 

=  D 

7.2 

= 

00 

• 

1 — 1 

P 

3.0 

-0.9 

5)  -  0.9  ,  p2  =  -0.9  ,  P3  =  -0.9  , 

P4  =  0 . 7  ,  P5  ~  -0 . 3  ,  Pg  ~  1.3. 

P5  and  p6 


In  determining 


» 


we  have  assumed  that  e  is 
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insignificantly  small  compared  to  0.3.  We  have 
the  equality  sign  by  the  symbol  "~M  whenever 
process  has  been  carried  out. 


6 )  | pr |  =  max  | p . |  j  a  =  6  . 

l<i<6  1 


7)  |  p6 1  >  p  .  6*1. 


8)  u  =  1  . 


9)  CX1,X2>X3]  =  [1 ,2+e , 1]  D  -  [i(l+e),  -|(l+e) 


10) 


max  u 
l<i<3 


k  =  3  . 


11)  D 


1 

9+5e 


1+e  3+e  -2 

3-3  5 


replaced 

similar 


9+5£  )  ]  . 


-5-3e  3+2e 


1 
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12)  I  =  [1,2,6]  . 


Step  2 


4) 


5) 


6) 


7) 


X1 

1 

10 ( 9+5e ) 

1^5  + 

1  CJ 

1 — 1 
-^r 

i — 1 

x2 

272 

P 

-59  - 

63e 

.59  .  . 

59 

P1 

90  s  p2 

90  5 

p3 

p4  ' 

217 

379 

p6  ' 

90  J  p5 

90  J 

1  Ps 

=  max  Ip. 
l<i  <6 

;  Ot  = 

5  . 

K 

>  P  .  5^1 

• 

147 

59 

90 


8)  u  =  -1  . 
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9)  [A 


u 

10)  - 

For  e 
ratios 

11)  D 


]_  j  A2  ,  A^  ] 


=  [2+e,4,-l]  D 


[15+6e+e2, 


-9+3e+e  ,  15-2e]  . 


X1  =  15+6e+e2  t  u  X2  =  9-3e-e2  .  u  X3 


^3 1  5+3e  s  d^2  3+2e  ’  d33 


=  -15+2e  . 


>  0  and  sufficiently  small,  the  maximum  of  the  above 
u  A, 


is 


-3 —  .  Therefore,  k  =  1 
d31 


1 


(9+5e) (15+6e+e  ) 


(l+e)(9+5e)  54+39e+5e  -45-25e 

2(9+5e) 


-27-2  4e-5e2  45+34e+5e2 


(_5_3e) (9+5e)  36e+29e2+5e2  90+4le-5e2 


12)  I  =  [5,2,6]  . 


69 


Step  3 


4) 

xi 

1737  +  26l2e  +  915e 

xp 

l 

p 

3204  +  1762e  -  10e2 

10(9+5e)(15+6e+e^) 

p 

1305  -  310e  -  7e2 

2637 

,  1305 

891 

1350  » 

p2  1350  » 

p3  ~ 

1350 

2817 

n  1305 

p6  ~ 

1305 

1350  5 

p5  ~  "1350  5 

1350 

6 )  |  p  J  =  max  |  p,  |  ;  a  =  4  . 

l<i<6  -1 


7)  |p4|  >  P  .  4^1. 


360e 


8)  u  =  1  . 
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9)  [X1,X2jA3]  =  [1,-1, 1]  D 

1 

=  - p~  x 

( 9+5e ) ( 15+6e+e  ) 

[(-6-2e) (9+5e) ,  8l+99e+39e2+5e3 ,-l8e-10e2 ] . 


10) 


u  Xp 
d32 


X. 

max  u  -5 — — 
l<i<3  31 


k  =  2  . 


11)  D  = 


1 

3+e 


0 

1 

3  +  e 


6+e 

(3+e)2 

1 

"3  +  e 


4e+e2 

(3+e)2 


(3+e)2 

1 

3+e 
6  +  8 

(3+e)2 


12)  X  =  [5,^,6]  . 
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Step  4 


4) 

xi 

= 

l83+121e 

10(3+e)2 

x2 

6 

3+e 

87-e+10e2 

P 

10 ( 3+e ) 2 

5) 


Pi  -  — 


75 

81 

n  93 

90  1 

p2  ~  "90  > 

p3  ~  "90  * 

87 

n  87 

87 

90  » 

p5  ~  "90  s 

p6  ~  90  ‘ 

6 )  I  p  ~  I  =  max  I p .  I  ;  a  =  3  . 
8  l<i<6  1 


7)  |p^|>p.  3  £  I  . 


8)  u  =  -1  , 


) 


6 


3— 4e-e^ 

(3+e)2 


10)  2 


u  A1 


max  u  -? — 
l<i<3  3i 


k 


11)  D  = 


1 

2 


0 


1 

3+e 


1  1+e 

2  6  +  2g 


0 

1 

3+e 

1 

3  +  e 


12)  I  =  [3,^,6]  . 

Step  6 


4) 

xi 

= 

2 

x2 

6 

3+e 

P 

6-2e 

6+2e 
-  - - 

*  ? 

( 3+e ) 


=  1  . 
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5)  px  -  0.9  , 

P2  ~  -0.8  , 

0 

. 

i — 1 

1 

1 

on 

Q. 

O 

• 

1 - 1 

^r 

Q- 

~  -0 . 9  , 

• 

0 

• 

f — 1 

? 

vo 

Cl 

6)  |pQ|  =  |  pi,  |  =  |  pr  I  =  max  Ip.  I  ;  a  =  3  or  4  or  6  . 

^  l<i<6  1 

7 )  I I  =  I P4 I  -  | Pg I  >  P  >  a  e  I  . 

The  last  line  shows  that  the  critical  subsystem  has  been 

reached.  To  find  the  minimax  solution  of  (4.1)  we  let 
e  -*■  0  to  obtain  x^  =  x^  =  2  and  p  =  1  . 

The  technique  shown  above  is  not  convenient  for  use 
on  a  computer  since  e  has  been  left  as  a  variable  in  the 
computations.  The  same  method  employed  on  a  computer  would 
require  the  use  of  symbol  manipulation,  which  is  not 
readily  available.  The  only  recourse  would  seem  to  be  to 
introduce  a  sequence  of  e’s  say  ei  >  e2  >  •••  en  ^ 
and  examine  what  happens.  If  no  significant  changes  are 
observed  after  a  certain  ,  then  one  could  accept  the 

minimax  solution  to  the  problem  with  appearing  as  a 

good  approximation  to  the  minimax  solution  of  the  original 
prob lem . 
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The  problem  of  finding  the  minimax  solution  of  (4.1)  was 
attacked  using  the  computer  program  MINMAX  given  in  the 
Appendix.  This  program  assumes  that  Haar’s  condition  is 
satisfied.  The  correct  solution  was  obtained  using  both 
single  and  double  precision  arithmetic,  using  as  initial 
guess  I  =  [1,2,3]  .  The  computations  did  not  break  down 
(as  in  the  exact  calculation  shown  in  Section  4.1)  because 
roundoff  errors  introduced  during  computations  yielded  very 
small  numbers  rather  than  zeros  at  the  critical  stages. 

In  effect,  the  roundoff  errors  perturbed  the  original 
system,  which  did  not  satisfy  Haar’s  condition,  into  an 
adjacent  system  that  did.  We  present  some  of  the  results 
of  the  single  precision  computations  to  3  significant 
digits.  (The  initial  subsystem  is  given  by  I  =  [1,2,3]  ). 

Step  2 

1)  T±  =  -2,  T2  =  1,  T3  =  1  . 

2)  s  — 1 ,  s 2  1,  s^  1  • 


' 


3) 


~i 


o 

• 

o 

0.5 

-0.5 

0.5 

-0.75 

1.25 

0.5 

o 

ro 

U1 

0.25 

4) 

X1 

2.10 

X2 

= 

o 

oo 

• 

1 — 1 

P 

-0.90 

0.9  , 

o\ 

• 

o 

1 

II 

C\J 

a. 

ii 

oo 

Q. 

-0.9 

0.7  , 

P5  =  -0.3  , 

P6  = 

1.30 

6)  | pg |  =  max^  | pi |  ;  a  =  6  . 


7)  | P6 1  >  P  .  6  t  I  . 


8 )  u  =  1  . 
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9)  CX1,A2,A3]  =  [0.5,-0.75,2.25]  . 


10) 

9 . 0  = 

u  A~ 

A. 

max  u  -5 —  ;  k 
l<i<3  31 

11) 

0.111 

0.333 

-0 . 222 

D  = 

0.222 

-0.333 

0.556 

-0.556 

0.333 

0  .  Ill 

k  =  : 


12)  I  =  [1,2,6]  . 


Step  2 


A) 

X1 

1 — 1 

VO 

0 

1 — 1 

1 _ 

x2 

= 

CM 

O 

OO 

P 

-0 . 656 

5)  p2  =  0.656  ,  p2  =  -O.656  ,  p^  = 
P4  =  2. 41  ,  p^  =  -4.21  ,  p6  = 


-1.63  , 

-0.656  . 
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6 )  | p  |  =  max  | p . |  ;  a  =  5  . 

D  l<i<6 


7)  |  p ^  |  >  p  5  £  I  . 


8)  u  =  -1  . 


9)  [X1,X2,X3]  =  [1. 67,-1.0,1.67]  . 


u  X  p  A . 

10)  3*0  =  — -t —  =  max  u  -? —  ;  k  =  2  . 

a32  l<i <3  a3i 


11) 


D 


0 . 667 

-0.333 

0.333 

-0.333 

0.333 

0 . 212x10 

-0 .212*10-6 

-0.333 

0.667 

12)  I  =  [1,5,6]  . 


' 
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Note  in  stage  (11)  that  ^  0  .  In  the  exact  computation, 

d^  =  0  .  The  computations  therefore  do  not  break  down. 

The  algorithm  then  considers  successively  the  subsystems 
[2>5,6],  [3,5,6],  and  [3,4,6],  It  reports  that  the  last 
subsystem  is  the  critical  one  and  gives  the  minimax  solution 
as  x^  =  x2  =  2.0  and  the  associated  minimax  error  as 

p  =  1 . 0  . 

This  result  is  somewhat  encouraging.  It  should  be 
noted,  however,  that  the  reverse  process  is  also  possible: 
that  is,  it  is  possible  for  roundoff  errors  to  perturb  a 
system  satisfying  Haar’s  condition  into  an  adjacent  one 
that  violates  Haar’s  condition.  Individual  cases  will 
probably  require  individual  treatment. 


. 


CHAPTER  V 


MINIMAX  APPROXIMATIONS  ON  A  CONTINUUM 


5 . 1  Introduction 

In  the  previous  two  chapters,  we  solved  the  problem 
of  determining  a  vector  c  =  (c^,...,c  )  to  minimize 


n 


M-,  =  max  |  f  (x.  ,y  .  )  -  l 
l<i <m  1  1  j=l 


I  °jSJ(xi.yi)l  » 


(5.1) 


where  (x^jy^)  ,  i=l,...,m  ,  were  distinct  points  on  the 

contour  C  .  Intuitively,  we  feel  that  if  the  number  of 
points  becomes  large,  leaving  no  wide  gaps  along  C  ,  the 
minimax  error  associated  with  minimizing  (5.1)  should  be 
a  good  approximation  to  the  minimax  error  associated  with 
minimizing 


n 


(5.2) 


Mp  =  max  |  f  (x  ,y  ) 
(x,y)eC 


I  c.g  (x,y) 
j=l  J  J 


Thus,  one  approach  to  finding  the  minimax  error  on  the 
curve  C  would  be  as  follows:  Compute  the  minimax  error 
for  m  points.  Double  the  number  of  points  and  again 


V  AMPUJISj  '$  3  l: 
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calculate  the  minimax  error.  If  the  minimax  error  does 
not  change  significantly,  then  stop.  Otherwise,  once  more 
double  the  number  of  points  and  continue  as  before. 

A  more  elegant  approach  was  first  conceived  by  Remes 
in  1934.  The  statement  of  the  algorithm  and  the  proof 
of  its  convergence  is  from  Cheney  ([3],  pp.  95-96). 

5 . 2  Remes  Algorithm 

This  algorithm  only  requires  that  the  functions 

be  continuous  on  the  curve  C  .  We  wish  to 
determine  a  vector  c  =  (c^,...,c  )  for  which  the  deviation 


n 


I  c.g  (s) 
0=1  J  0 


(5.3) 


M( c )  =  max  | f ( s ) 
seC 


is  a  minimum.  We  may  assume  without  loss  of  generality 
that  gia...jgn  are  linearly  independent,  for  if  they  are 
not,  we  can  replace  them  by  a  smaller  number  of  functions 
that  are  independent  without  raising  the  minimum  deviation 


(5.4) 


p  =  min  M(c)  . 
c 


We  define  a  residual  function  r(c,s)  by 
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(5.5) 


n 


r(c,s)  =  f(s)  -  l  cg.(s) 

j=l  J  J 


■f”  K 

Then  Remes  algorithm  is  as  follows:  At  the  k  step  we  are 

!/■ 

given  a  finite  subset  S  of  C  .  Select  a  coefficient 
vector  c  to  minimize  the  function 


(5.6) 


M  ( c )  =  max  | r ( c  ,s  ) 

Qk 

seS 


This  can  be  accomplished  by  the  procedure  given  in  Chapter 

k  i  1c  i 

III.  Choose  s  eC  to  maximize  |r(c  ,s)|  .  Thus, 


(5.7) 


r(ck,sk) |  =  M(ck)  . 


k+ 1  V  \c 

Now  start  again  with  the  finite  set  S  =  S  u{s  }  .  At 
the  beginning,  S1  may  be  arbitrary  except  that  the 
set  of  n-tuples  §  =  [g^(s ) , . . . ,g  (s ) ]  corresponding  to 
seS1  should  be  of  rank  n  .  The  algorithm  stops  when 

|r(c  ,s  ) |  is  insignificantly  different  from  min  max 

c  0k 

seS 

|r(c,s)|  .  A  theorem  associated  with  Remes  algorithm  is 


k  k 

Theorem  5.1.  Urn  M  (a  )  =  p  .  The  sequence 

k-H=° 

k 

c  is  hounded  and  its  cluster  points  minimize  M  . 


. 


.  , 


J92 
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Proof:  We  define 


(5.8)  II  c||  =  l  |  c  |  . 

j=l  J 

It  can  be  shown  that  this  function  defines  a  norm  in  the 
linear  space  .  By  our  assumption  on  S'*"  ,  it  follows 

that 


(5.9) 


<f> 


min 

II  e|  1=1 


max 

seS"*- 


n 

I  c.g  (s) |  >  0  . 
j=l  J  0 


Consequently , 


M1(c) 


max  |f(s) 
seS"*" 


l  c  g  (s)| 
j  =  l  J  J 


=  II  f  -  I  c.g  || 
j=l  J  J 


(5.10) 


*  II  l  C  4  S4  II 
j=l  J  J 


=  max 
s  eS^ 


l  c  g  (s)| 
j=l  0  J 


max  | f ( s ) |  . 

seS"*" 


' 
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Now  let 


(5.11) 


K  . 


Then 


(5.12) 


max 
1 


seS 


n 

I  l  c  g  (s) |  =  K 
j=l  J  J 


max 
1 


seS 


n 

I 

j=i 


Cj/K 


gj  t s )  i 


Therefore,  we  may  rewrite  (5.10)  as 


(5.13)  M  (c)  >  ||c||,  max  |  £  c.g  (s)|  -  max  |f(s)|  , 

x  Q1  j=l  J  J  C1 

seS  d  seS 


and  hence  by  (5.9), 


(5.1^)  M1(c)  >  ||  c ||  <}>  -  max  |f(s)|  . 

s  eS1 


Now  suppose  that  ||c||,  >  2  max  |  f* ( s )  |  /<(>  .  Then 

seS1 
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Mk(c)  >  M1(c) 


(5.15) 


>  max  | f (s ) 
seS1 


>  Mk ( 0 )  , 


where  0  is  the  zero  vector.  Therefore,  c  does  not 

enter  the  competition  to  minimize  any  of  the  functions 
k  k 

M  .  Thus  the  sequence  c  generated  by  the  algorithm 
is  bounded.  Now  we  have 


(5.16) 


c  S 


k+1 

C 


C 


> 


and  hence 


(5.17)  Mk(c)  <  Mk+1(c)  <  M(c) 


for  any  vector  c  .  Therefore, 


..k,  kx  ^  ..k+1,  k+lx  ^ 

M  (c  )  <  M  (c  )  <  p 


(5.18) 
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Thus  for  some  e  ^  0  ,  M^(c^)  -+■  p  -  e  .  We  must  show  that 
e  =  0  .  Since 

(5.19)  | r(b  ,s )  -  r(c,s)|  =  |  l  (c  -b.)  g,(s)| 

j=l  J  J  J 


*  T  l|t>-c||  1  , 


where 


(5.20)  T  =  max  max  |g.(s)|  , 

l<j <n  seC  J 


it  follows  that 


(5.21)  |r(b,s)|  <  |r(c,s)|  +  T  (b-c^ 


and  furthermore 


(5.22) 


M(b)  <  M(c)  +  T  ||  b— c  ||  1  . 


Suppose  now  that  e  >  0  .  Let  b  denote  any  cluster  point 


(  '  .  •  - 
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\r 

of  the  sequence  c  .  For  any  6  >  0  ,  we  may  find  an 
index  k  such  that 


(5.23) 


||b-ck||1  <  6 


and  an  index  i  >  k  such  that 


(5.24)  || b-c^||  ^  <  6  . 

Then 

(5.25)  ||o1-clc||1  £  26 
and 

p  <  M(b)  <  M(ck)  +  T6 

=  | r(ck,sk) |  +  T5 

(5.26)  <  |r(ci,sk)|  +  3T6 

<  Mi(ci)  +  3T6 

<  p  -  e  +  3T6  . 


.  a  -  q  ■  2 
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Picking  6  small  enough  so  that  3T6  <  e  yields  a 
contradiction.  Therefore,  e  =  0  and  (5.26)  can  be 
rewritten 


(5.27)  p  <  M(b )  <  p  +  3T6  . 

Letting  6+0  proves  the  assertion. 

Remes  algorithm  requires  that  we  locate  the  extremum 
of  the  residual  function  r(c,s)  .  This  can  be  a  difficult 
problem  even  in  well-behaved  cases.  There  would  be  no 
problem  in  constructing  an  example  to  defeat  any  particular 
method . 

We  have  chosen  the  obvious  and  simple  method  of 
calculating  the  residual  on  a  discrete  set  of  points  chosen 
in  advance.  Originally,  this  discrete  set  was  chosen  as 
a  set  of  equally  spaced  points.  Examination  of  the  residual 
function,  however,  has  shown  that  it  has  approximately  the 
form  of  a  Chebyshev  polynomial  ([13],  p.  60),  with  the 
extrema  crowded  towards  the  ends  of  the  range  under 
consideration.  Therefore,  the  discrete  point  set  was  chosen 
to  be  the  zeros  of  a  Chebyshev  polynomial.  As  a  rule  of 
thumb,  if  we  are  approximating  f  by  n  functions 
gl,*’*,gn  5  then  the  discrete  point  set  was  chosen  to  be 
the  zeros  of  the  (20n)th  Chebyshev  polynomial,  although 
no  significant  change  in  the  minimax  error  was  evident  if 


j  :)  .  ■  >: 


( cf)M  fe  q 


.  (a  o j*i 


. 
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we  took  the  discrete  point  set  to  be  the  zeros  of  the 
Chebyshev  polynomial  of  (10n)th  degree  or  higher. 


CHAPTER  VI 


NUMERICAL  RESULTS 

The  present  experiments  were  performed  using  the  IBM 
System  3 6 0  FORTRAN  IV  source  language  and  the  IBM  360/67 
computer  system  at  the  University  of  Alberta,  Edmonton. 

6 . 1  Test  Problems  and  Results 

For  clarity,  we  shall  elaborate  the  meaning  of  several 
terms  which  occur  in  the  tabulation  of  the  data. 

Boundary  C:  The  boundary  of  the  region  of  space  over  which 

the  solution  is  sought. 

Boundary  values :  A  continuous  function  f  defined  on  the 

boundary.  If  there  are  more  than  one,  then  the 
functions  will  be  labelled  (a),  (b),  ...  . 

Approximating  functions:  Harmonic  functions  selected  from 

a  set  of  harmonic  functions  which  is  complete 
over  the  space  of  functions  which  are  continuous 
on  the  boundary.  Usually,  these  functions 
will  be  selected  from  the  set  of  harmonic 
polynomials  (see  Section  2.3). 

Boundary- value  problem:  This  will  be  either  the  Dirichlet 

problem  or  the  mixed  problem  for  Laplace’s 
equation  (see  [16],  p.  169,  for  the  definition 
of  the  mixed  problem) . 
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Minimax  error,  :  The  computed  minimax  error  when  n 

approximating  functions  are  employed. 

Minimax  solution:  This  will  be  written  in  the  form 

u  ~  cpS]_  +  .  .  .  +  cngn  »  where  u  is  the 
required  solution  of  the  Dirichlet  problem, 
gl,...,gn  are  the  approximating  functions, 
and  c1,...,cn  are  the  computed  coefficients 
that  minimize  the  maximum  error  on  the  boundary. 
Least-Squares  solution:  This  will  be  of  the  same  form  as 

the  minimax  solution  except  that  the  coefficients 
minimize  the  sum  of  the  squares  of  the  errors 
at  discrete  points  along  the  boundary. 

Error  at  isolated  interior  point(s):  The  absolute  value  of 

the  difference  between  the  theoretical  value 
and  the  computed  value  at  the  given  point (s). 

This  can  be  given  only  if  the  true  solution 
is  known. 

6.1.1  Problem  1  For  complete  details  see  Davis  [12], 
Boundary ;  A  bean-shaped  region  was  obtained  from  a  free¬ 
hand  drawing  on  coordinate  paper  (see  Figure 
6.1).  The  boundary  is  ’’defined"  by  means  of 
43  points  (see  Table  6.1).  These  points  are 
not  distributed  equally  on  the  boundary;  more 
points  appear  where  the  curvature  is  greater. 


' 
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Figure  6.1  Bean-shaped  region  of  Problem  1 
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Table  6.1 


Points 


defining  the  bean-shaped  region 


Pt »  no. 


Abscissa  Ordinate 


1 

.  000 

.  110 

2 

-.050 

.  108 

3 

-.100 

.115 

-.160 

.150 

5 

-.220 

.205 

6 

-.320 

.300 

7 

-.400 

.358 

8 

-.500 

.  420 

9 

-.550 

.436 

10 

-.600 

0 

on 

• 

11 

-.644 

.  400 

12 

- »  660 

.350 

13 

-.655 

.300 

14 

-.635 

.  200 

15 

-.595 

.  100 

16 

-.552 

.  000 

17 

-.500 

-.105 

18 

-.440 

-.200 

19 

-.400 

-.250 

20 

-.350 

-.300 

21 

-.300 

-.344 

22 

-.204 

-.400 

23 

-.100 

-.436 

24 

.  000 

-.448 

25 

.  100 

-.442 
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Pt 


Table  6.1  (Continued) 


.  no . 

Abscissa 

Ordinat 

26 

.230 

-.400 

27 

.300 

-.350 

28 

.353 

-.300 

29 

.430 

-.200 

30 

.477 

-.100 

31 

.510 

.  000 

32 

.522 

.100 

33 

.520 

.  160 

34 

.500 

.240 

35 

.456 

.300 

36 

.  400 

e  3  3  0 

37 

.360 

.337 

38 

.300 

.320 

39 

.250 

.290 

40 

.200 

.245 

41 

.  150 

.  200 

4  2 

.100 

.  160 

43 

.050 

.128 

Boundary  values: 


(a)  f 

(b)  f 

(c)  f 

(d)  f 

(e)  f 

(f)  f 

(g)  f 

(h)  f 


0  a  x  <  0 ;  f  =  x,  x  >  0  , 


0  j  x  <  0 ;  f  =  x  j  x  >  0 


0  a  x  <  0 ;  f  =  x  a  x  >  0 


0  a  x  <  0 ;  f  =  x  a  x  >  0 


0  a  x  <  0 ;  f  =  x  a  x  >  0  . 


2  2 

x  +  y  (Torsion  problem). 


2  2 
exp ( x  -  xy  +  2y  ) 


ln[ ( x  +  1. 5)2  +  y2]  . 


(i)  f  =  ln[x  +  y  +  (y-1)2]  . 

(j)  f  =  ex  cos  y  +  ln[x2  +  (y-1)2]  . 

n  n 

Approximating  functions:  1,  Re  z  ,  Im  z  ;  n=l,2. 
Boundary  value  problem:  The  Dirichlet  problem. 
Minimax  error:  See  Table  6.2. 
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Minimax  solution: 

This  is  for 

function 

(j)  for  n 

=  5  . 

u  ~  1.00244 

+ 

.997007 

Re 

z 

— 

1.99034 

Im  z 

+ 

1.47952 

Re 

2 

z 

- 

. 0110482 

Im  z2 

+ 

.187437 

Re 

z3 

+ 

.628927 

Im  7? 

- 

. 3^5989 

Re 

4 

z 

+ 

.0152750 

Im  z 4 

- 

d33709 

Re 

z5 

- 

.322899 

Im  z5 

Least  squares  solution:  The  coefficients 

are  quoted 

from 

[12] 

• 

This  is  for  the  : 

same 

problem  as 

above . 

u  ~  1.00173 

+ 

.997339 

Re 

z 

- 

1.99119 

Im  z 

+ 

1.48065 

Re 

2 

z 

- 

.  00949996 

Im  z2 

+ 

.188957 

Re 

z3 

+ 

.623677 

Im  z^ 

— 

.355600 

Re 

4 

z 

+ 

,  024526 

,  4 
Im  z 

_ 

. 11960 

Re 

z5 

— m 

.28034 

Im  z5 

For  this  approximation  the  maximum  error  on 
the  boundary  is  .0069. 
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Table  6.2  Results  for  bean-shaped  region 


Minimax  error  for  function 

n 

(a) 

(b) 

(c) 

(d) 

(e) 

1 

.  142 

.0823 

.0477 

.  0266 

.0145 

2 

.0870 

.0410 

.  0232 

.0135 

.00776 

3 

.0725 

.  0284 

.0131 

.00657 

.00357 

4 

.0613 

.0234 

.  00961 

.00378 

.  00156 

5 

.0562 

.0182 

.00757 

.00321 

.00138 

6 

.0480 

.0164 

.00662 

.00276 

.00123 

7 

.  0446 

.0141 

.00546 

. 00233 

.00105 

8 

.0414 

.0126 

. 00483 

. 00202 

.000909 

9 

.  0382 

.0113 

.00416 

1 — 1 

O 

O 

.000783 

Minimax  error  pn  for  function 

n 

(f) 

(g) 

(h) 

(i) 

(J) 

1 

.190 

.544 

.613 

.384 

.277 

2 

.148 

.338 

.  0322 

.192 

.0858 

3 

.  0878 

.184 

.00930 

.  0940 

.0258 

4 

.  0707 

.116 

.00250 

.0429 

.0110 

5 

.  0599 

.0964 

.000759 

.  0216 

.  00434 

6 

.0505 

.  0810 

. 000226 

.0122 

.00179 

7 

.  0433 

.  0698 

.0000707 

.00879 

.000738 

8 

.0371 

.  0607 

.0000204 

.00645 

.000321 

9 

.0298 

.0487 

. 000000545 

.00491 

.000134 
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Error  at  isolated  interior  points:  The  results  (see  Table 

6.3)  pertain  to  function  ( j )  for  both  the 
minimax  and  least-squares  solution. 


Table  6.3  Discrepancies  at  points  interior 

to  bean  along  y  =  0 


X 

minimax  discrepancies 

least-squares 

discrepancies 

-.5 

o  0033 

.  0014 

-.4 

.  0021 

.  0010 

-.3 

.  0020 

.  0011 

-.2 

.0023 

.  0015 

-.1 

.  0025 

.  0017 

0 

.  0024 

.0017 

.1 

.  0020 

.0013 

.2 

.0013 

.  0007 

.3 

.  0006 

.  0001 

.4 

.  0005 

.  0009 

6.1.2  Problem  2 
Boundary :  The  ellipse 


2  ...  2 
x  +  4y 


1  . 
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Boundary  values : 

(a)  f  =  03  x  <  0;  f  =  x,  x  >  0  . 

(b)  f  =  0,  x  <  0;  f  =  x2,  x  >  0  . 

(c)  f  =  0,  x  <  0;  f  =  x^ ,  x  >  0  . 

(d)  f  =  ln[(x  +  1.5)2  +  y2]  . 

Approximating  functions:  1,  Re  zn;  n=l,2a...  .  Only  the 

real  parts  of  zn  contribute  since  the  boundary 
and  boundary  functions  are  symmetric  with 
respect  to  the  x  axis.  Thus  only  the  top 
half  of  the  boundary  need  be  considered. 

Boundary-value  problem:  The  Dirichlet  problem. 

Minimax  error:  See  Table  6.4. 
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Table  6.4  Results  for  the  ellipse 


Minimax 

error  p 

Kn 

for  function 

n 

(a) 

(b) 

(c) 

(d) 

1 

.250 

.281 

.317 

.385 

2 

.0577 

.0856 

.134 

.139 

3 

.0525 

.0138 

.0370 

.0577 

4 

.0311 

.0138 

.00442 

.0258 

5 

.0308 

.0531 

.  00444 

.0120 

6 

.  0210 

.0531 

.00133 

.00571 

7 

.0209 

.00277 

.00133 

.00278 

8 

.0159 

.  00277 

,000566 

.00137 

Exact  value  of  solution  for  function  (d)  at  the  point 

(0,0)  :  u( 0 , 0  )  =  In  2.25  =  0.8109  . 


Computed  value  of  approximating  solution  at  the  point 

(0,0)  for  n  approximating  functions:  See  Table  6,5. 
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Table  6.5  Computed  solution  for  function 

(d)  at  (0,0) 


n 

computed  value 
at  (0,0) 

discrepancy  between  computed 
and  true  solution 

2 

.6105 

.  2004 

3 

.  8544 

.0435 

4 

.  8207 

.  0098 

5 

.  8077 

.  0032 

6 

.8103 

.  0006 

7 

.  8112 

.0003 

8 

.  8110 

.  0001 

9 

.8109 

.  0000 

6.1.3  Problem  3 


Boundary :  Square  given  by  the  lines  x=0,  x=l,y=0, 

y  =  1  . 


Boundary  values :  (a) 


1,  y  =  I,  0  <  x  <  1 


f  ■  s 


0,  otherwise 


(b) 


x  -  x ,  y  =  1 ,  0  <  x  <  1  . 


=  < 


0,  otherwise  . 


, 
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Approximating  functions:  sinh  irny  sin  imx,  n=l,2J...  . 

These  functions  are  harmonic  and  satisfy  the 
homogeneous  boundary  conditions.  Because  of 
the  symmetry  of  the  problem  along  the  line 
x  =  i  ,  the  approximating  functions  used  are 
sinh  7r(2n-l)y  sin  ir(2n-l)x  ,  n=l,2,...  . 

Minimax  error:  See  Table  6.6. 


Table  6.6  Results  for  Square 


Minimax  error  p 

Kn 

for  function 

No.  of  functions 

(a) 

(b) 

2 

1.00 

.0024 

3 

1.00 

.  0010 

4 

1.00 

.00054 

5 

1.00 

.  00034 

6 

1.00 

.  00023 

7 

1.00 

. 00016 

8 

1.00 

. 00012 

9 

1.00 

.  00010 

Exact  solution  for  function  (a)  at  (0.5, 0.5):  u(0.5,0.5) 


0.25  (see  [22],  p.  105). 


Computed  value  of  solution  at  (0.5, 0.5)  for  n 

approximating  functions:  See  Table  6.7. 


102 


Table  6.7  Computed  solution  for  function 

(a)  at  (0. 5,0.5) 


n 

value 

2 

.2474 

3 

.2519 

4 

.2505 

5 

.2504 

6 

.2504 

7 

.'250  4 

8 

.2455 

9 

.2501 

6.1.4  Problem  4  For 
Boundary :  Rectangle  given 


complete  details, 
by  the  lines  x  = 


see  Charmonman  [2]. 
0,  x  =  -2, 


y  =  0,  y  =  -1  . 


Boundary  values : 


x,  y  =  -1,  -1.5  <•  x  <  0  . 


f  =  < 


0 

V 


> 


otherwise  . 
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Boundary-value  problem:  The  mixed  problem.  The  normal 

derivative  of  the  solution  is  specified  along 
y  =  -1,  -2  <  x  <  -1.5  .  Elsewhere,  the  value 
of  the  solution  is  given. 


Approximating  functions:  i)  sinh  sin  ,  n=l,2,... 

where  the  value  of  the 
solution  is  given. 


.  .  \  7m  ,  Tmy 
li )  — ^  cosh  —^2 


Tmx  ,  0 

sm  — ^  >  n=l ,  2  , . 


where  the  normal  derivative 


of  the  solution  is  specified. 


Minimax  error:  See  Table  6.8. 


Table  6.8  Results  for  Rectangle 


No.  functions 

Minimax  error  p 

2 

.7461 

3 

.6317 

4 

.  6l40 

5 

.  5875 

6 

.5264 

7 

.5204 

8 

.5052 

9 

.  4612 

10 

.4598 

. 
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6 . 2  Comments  on  Numerical  Results 

Problem  1:  Davis  [12]  mentions  this  problem  as 
a  difficult  test  case  because  of  the  nonconvexity  of  the 
domain.  Nevertheless,  we  were  interested  in  testing  our 
procedure  for  a  fairly  intricate  region. 

Boundary  functions  (a)  to  (e)  yield  a  family  of 
boundary  values  of  increasing  smoothness,  The  table  of 
minimax  errors  show  a  decrease  in  the  maximum  boundary 
error  as  the  boundary  functions  become  smoother.  In 
addition,  the  rate  of  decrease  of  the  minimax  error  as 
we  introduce  more  approximating  functions  is  greater  for 
smoother  boundary  functions. 

Function  ( j )  is  harmonic  over  the  region  over  which 
the  solution  is  required  and  hence  the  boundary  function 
is  the  solution  to  the  Dirichlet  problem.  This  test  case 
was  tried  in  [12]  using  a  least-squares  approach  mentioned 
in  Chapter  I.  We  note  that  the  maximum  error  along  the 
boundary  is  smaller  using  our  minimax  procedure  than  the 
maximum  error  given  by  Davis'  least-squares  technique. 
However,  the  discrepancies  in  the  interior  seem  to  be 
larger  for  the  minimax  technique.  In  both  approaches,  the 
discrepancies  in  the  interior  are  less  than  the  computed 
maximum  error  along  the  boundary. 


105 


Problem  2:  For  test  functions  we  selected  four 
of  the  functions  which  were  tried  in  Problem  1.  (Boundary 
functions  (a),  (b),  (c),  (d)  for  the  ellipse  correspond 
to  boundary  functions  (a),  (b),  (c),  (h)  for  the  bean.) 

For  the  family  of  functions  (a)  to  (c)  ,  we  see  that 
the  minimax  errors  for  the  ellipse  are  less  than  the 
corresponding  minimax  errors  for  the  bean.  This  is  probably 
due  to  the  nonconvexity  of  the  bean-shaped  region.  On 
the  other  hand,  the  minimax  errors  for  function  (d)  are 
bigger  for  the  ellipse.  This  is  a  result  of  the  singularity 
of  the  boundary  function  being  closer  to  the  ellipse  than 
to  the  bean.  For  function  (d)  ,  we  have  listed  the 
differences,  for  increasing  numbers  of  approximating 
functions,  between  the  true  solution  and  the  computed 
approximate  solution  at  an  interior  point.  The  differences 
become  smaller  as  the  number  of  approximating  functions 
increases.  We  note  also  that  the  discrepancies  are  well 
within  the  computed  maximum  error  along  the  boundary. 

Problem  3:  The  two  boundary  functions  used  show 
the  difference  in  results  between  a  case  where  the  boundary 
function  possesses  a  discontinuity  on  the  boundary 
(function  (a))  and  a  case  where  it  does  not  (function  (b)). 
The  solution  for  the  Dirichlet  problem  using  function  (a) 
can  be  written  as  an  infinite  series  whose  value  is  known 
at  the  interior  point  (0. 5,0.5)  .  We  have  shown  how  the 


, 


. ,  n  ;:oq  loi'xsrfni  as  jb  rtci^wl  .a  0 Jb  i  :o»  qq 
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computed  value  approaches  the  exact  value  at  this  point  as 
more  approximating  functions  are  introduced.  The  approach 
to  the  true  solution  is  quite  uniform  except  for  8 
approximating  functions  where  the  sudden  divergence  is 
something  of  a  mystery,  although  it  in  no  way  contradicts 
any  theory  that  we  have  presented. 

Problem  4;  Although  we  have  not  dealt  with  the 
mixed  problem  in  the  text  of  this  thesis,  we  were  interested 
to  see  how  our  minimax  procedure  would  fare.  For 
approximating  the  normal  derivative  of  the  solution  along 
the  boundary  we  have  used  a  standard  technique  in  numerical 
analysis:  If  L  designates  a  linear  operator  and  u  a 

function,  then  we  approximate  L(u)  by  L( approximation 
to  u  ) .  It  is  interesting  to  note  that  the  minimax  error 
still  decreases  as  the  number  of  approximating  functions 
increases.  This  leads  us  to  believe  that  the  technique 
may  be  used  with  success  in  other  problem  areas,  such  as 
the  Neumann  problem  for  Laplace’s  equation  and  problems 
involving  multiply-connected  domains.  The  large  minimax 
errors  are  probably  due  to  the  singularity  at  the  point 

(-1.5,-lV  • 
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CHAPTER  VII 


CONCLUSIONS  AND  SUGGESTIONS 
FOR  FURTHER  RESEARCH 


7 . 1  Conclusions 

The  following  conclusions  have  been  drawn  from  the 
numerical  computations: 

1)  As  opposed  to  finite-difference  techniques,  the  method 
presented  is  easy  to  apply  to  boundaries  of  irregular 
shape  (for  example,  the  bean-shaped  region).  Further¬ 
more,  obtaining  the  approximate  solution  at  any  desired 
point(s)  in  the  interior  of  the  region  is  relatively 
simp le . 

2)  The  technique  is  certainly  as  good  as  other  collocation 
techniques  such  as  the  point-matching  method  and  the 
least-squares  method.  The  maximum  error  is  usually  a 
strictly  decreasing  function  of  the  number  of  approxi¬ 
mating  harmonic  functions.  This  is  a  virtue  that  the 
point-matching  method  does  not  always  possess  (see 
Section  1.3.1).  The  maximum  error  along  the  boundary 
is  smaller  for  this  method  than  for  the  least-squares 
method,  and  hence  the  bound  of  the  error  between  the 
true  solution  and  the  computed  solution  is  smaller. 
Furthermore,  the  maximum  error  is  part  of  the  automatic 
output  of  the  routine,  thus  requiring  no  additional 
computation  to  obtain  it. 


■ 
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3)  The  problems  discussed  show  clearly  the  sensitivity 
of  the  method  to  both  the  boundary  functions  and  the 
geometry  of  the  boundary.  The  best  results  seem  to 
come  from  boundary  functions  which  are  regular  in  a 
large  portion  of  the  plane,,  Poorer  results  are 
obtained  for  boundary  data  of  low  continuity  and  for 
boundary  data  which  possess  singularities  near  or  on 
the  boundary.  The  method  prefers  convex  regions 
(such  as  the  ellipse)  over  nonconvex  regions  (such  as 
the  bean),  although  convexity  does  not  seem  to  play 
as  important  a  role  as  the  smoothness  of  the  boundary 
functions . 

7 . 2  Suggestions  for  Further  Research 

1)  The  heart  of  our  computational  algorithm  is  finding 

the  minimax  solution  of  an  overdetermined,  inconsistent 
system  of  linear  equations.  We  have  presented  only 
one  method  of  determining  this  minimax  solution, 
namely  the  exchange  algorithm,  Cheney  [3]  briefly 
discusses  two  other  methods  -  Polya's  algorithm  and 
the  Descent  algorithm.  A  paper*  has  come  to  our 


*  Bartels,  R.H.,  and  Golub,  G.H.,  "Stable  Numerical  Methods 
for  Obtaining  the  Chebyshev  Solution  to  an  Over¬ 
determined  System  of  Equations",  JACM,  11,  1968, 
pp.  401-406, 


o 
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attention  as  this  thesis  is  in  print  which  presents 
yet  another  technique.  The  method  uses  the  exchange 
algorithm  basically,  but  computes  the  inverse  anew 
each  iteration  rather  than  updating  the  inverse  from 
the  previous  iteration.  The  author  claims  the  method 
to  be  more  stable  than  the  exchange  algorithm  as  we 
have  presented  it.  Furthermore,  the  algorithm  requires 
only  that  the  rank  of  the  coefficient  matrix  be  n  , 
rather  than  that  the  rows  of  the  coefficient  matrix 
satisfy  Haar’s  condition.  A  comparative  analysis  of  these 
algorithms  should  be  attempted. 

2)  Certainly,  the  Dirichlet  problem  for  Laplace’s  equation, 
although  important  in  its  own  right,  is  only  one  problem 
in  the  theory  of  partial  differential  equations. 

Numerical  investigation  of  the  mixed  problem  and  the 
Neumann  problem  for  Laplace’s  equation  should  be 
attempted.  Other  areas  of  investigation  should  include 
elliptic  equations  (of  which  Laplace's  equation  is  a 
subset),  boundary-value  problems  which  involve  multiply- 
connected  domains,  and  the  biharmonic  equation. 

3)  We  have  given  no  attention  to  the  use  for  which  the 
solution  of  the  Dirichlet  problem  is  intended.  Many 
problems  require  the  integral  of  the  solution  over  a 
region.  For  these  problems,  it  might  be  better  to 


minimize 


110 


n 

M  =  /  |  f  (s  )  -  l  a.  g.  (s  )  |  ds 
C  j=l  J  J 


rather  than  the  maximum  error  along  the  boundary. 
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APPENDIX 


LISTINGS  OF  FORTRAN  IV  SUBPROGRAMS 

The  following  pages  contain  listings  of  the  IBM 
System  360  Fortran  IV  subprograms  used  in  our  numerical 
experiments . 

Subroutines  DECOMP,  SOLVE ,  IMPRUV,  and  INVERT  have 
been  adapted  from  a  set  of  programs  given  in  Forsythe  and 
Moler  [15].  Details  can  be  found  therein.  Nevertheless, 
for  completeness ,  we  give  a  brief  description  of  each 
program. 

Subroutine  DECOMP  (N , A, UL, DET , ISGND )  uses  Gaussian 
elimination  to  find  N-by-N  triangular  matrices  L  and  U 
so  that  LU  =  PA  ,  where  PA  is  the  matrix  A  with  its 
rows  interchanged <,  The  matrix  L  -  I  +  U  ,  where  I  is 
the  identity  matrix,  is  stored  in  UL  0  DET  is  the 
determinant  of  A  and  ISGND  is  the  sign  of  the  determin¬ 
ant  of  A  . 

Subroutine  SOLVE  (N,UL,B,X)  uses  the  LU  factorization 
from  DECOMP  to  find  an  approximate  solution  to  a  system 
of  equations  AX  =  B  . 

Subroutine  IMPRUV  (N, A, UL,B,X, DIGITS)  requires  a  copy 
of  the  original  matrix  A  ,  its  LU  decomposition,  a 
right-hand  side  B  ,  and  the  approximate  solution  X 
computed  by  SOLVE ,  It  carries  out  the  iterative  improve- 
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ment  process,  until,  if  possible,  X  is  nearly  accurate 
to  machine  precision  (about  6  digits  on  the  IBM  S/36O/67). 

Subroutine  INVERT  (N,A,AINV)  uses  DECOMP,  SOLVE,  and 
IMPRUV  to  calculate  the  inverse  of  A  ,  which  is  then 
stored  in  AINV. 

Subroutine  MINMAX  computes  the  minimax  solution  and 
the  associated  minimax  error  of  an  overdetermined, 
inconsistent  system  of  linear  equations  Ax  =  b  . 

Subroutine  REMES  computes  the  vector  c  =  ( c ^ , „ . . , c^ ) 
that  minimizes 


n 

M  =  max  | f ( x,y )  -  l  a.g.(x,y) 

(x,y)eC  j=l  J  J 


Subroutine  ERRMAX  finds  (x2y)eC  such  that 


n 


f(x,y)  -  I  c  g  (x,y 


j=l 


max  |f(x,y) 
(x,y ) eC 


n 

v 


L  c.g  (x,y) 
j  =  l  J  J 


n 

v 


Subprogram  U  computes  1  c.g.(x,y) 


j=l 


J  J 


For  the  sake  of  completeness,  we  have  included  the  sub¬ 


programs  SGN(X),  F ( X, Y ) ,  G ( K , X , Y ) ,  and  BOUND(X)  .  SGN(X) 
is  the  program  for  the  mathematical  function  sgn  x  „ 
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F(X,Y)  is  the  program  for  the  function  f  specified  along 
the  boundary;  G(K,X,Y)  is  the  program  for  the  approxi¬ 
mating  function  g^  ;  BOUND(X)  is  the  program  of  the 
boundary  under  consideration.  These  last  three  programs 
are  valid  only  for  problem  2,  function  (d)  . 

Subroutine  SING  is  used  by  other  routines  to  indicate 
the  occurrence  of  an  error  condition. 


m  ' 
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o  UljkCU  T  i  .4  c  iJ  £  CCriP  (  i\ »  A  t  ul  f  J  c  1  ,l3GNUi 

L  I  Fi  L  i\  3  I  L  iN  A  1  2 0  »  2 0  I  >  C  l  (  2  0 , 2  0  J  »  A L  a t_  t  A  {  2  0  J  ,  I  2  A  i  20  i 

CUMKuN  I  HA 

G  U  3  I  =  i  ,  N 

i  P  S  (  I J  =  I 

R  C  k  !\  R  *M  =  C  .C 

L  C  2  J  —  1  ,  hi 

llL  t I »  J  J  =  A  1 i ,  J ) 

2  k  C  rt  N  K  M = A  rt  m  A  i  l  K  U  u  t\  k  M  i  A  b  d  (  U  L  l  I  ,  J  )  )  J 
IF  (  k Li  w f\ K M  i  d  i  ^  i  j 

3  SCAL£SiIJ=i.O/ruj'wNKM 
Gb  I C  3 

4  CALL  3  Ikb  l  1 ) 

Kb  TURN 

o  CuNT INuE 
0  £  I  =  1 .  C 
IS-1 
NM1=N-1 
OU  io 
b  I  b  =  C  •  C 
Cu  ii  I  =  K»i\ 

ip=ipa  i  n 

blbb^AoolCLi  lP,K)#3CALr:3l  I  P  J  ) 

It-(Ai^L-dioJ  1  1  i  1  i  ,  10 
1G  biG  =  3IZc: 

I G  X  P  IV  —  I 
il  CONTINUE 

IF l bio ) 13, 12 , 13 
Il  lAlL  3 1 N  G  (  2  J 
Rt  TCKiN 

13  IF  1  IGXPI V~Ki  14, 13 , 14 

14  IS-- IS 
J=  I  V  3  I  k  j 

IP  S ( K  J  =  i P  3 ( IlXPIV) 

IPbl  I  u  X  P  l  V  I  —  J 
Id  KP= I P  3 ( K ) 

PI VCJ=CL ( RP, KJ 
DtT=CbT*PI VCT 
Rp  1  =  in+  1 
CU  lo  I=KP1»N 
IP= IPS  (  I  ) 

EM=-UL  l  i  P  ,  R  ) /PIVUT 
CL  l  1 P  i  in  )  —  ~LFi 
CU  16  J-  is  P  1  ,  h 

lo  G  L  (  I  P  ,  J  )  =  G  L  I  IP,J)  +  tP>^UL(KP»J) 

K  P  —  I  P  3  {  N  ) 

DtI~i35i;CbT*Ui.(KP»N} 

IGGNC^SGinIulT  ) 

I F  (  b  L  t  K  P  ,  N  )  )  1 9  ,  i  d  ,  L  V 
Id  CALL  3  i  No ( 2  J 


19  HE 1  URN 
tl\u 

' 

. 


o  U  b  K  G  U  1  1 1\  £  SoLVb((\tOL»b»X) 

DlfocNbiUN  ULI2G»20)  j  t_>  (  2  0  )  »  X  i 20)  » I P S ( 2 0 ) 

C  G  jM  ivi  u  N  IPG 

l\P  1  =  N+  i 

lp-ipsm 

XU)=BUP) 

C  U  2  I  =  2  T  N 
IP-IPb (  i  ) 

IM  1=1-1 

surc-c.c 

DO  1  J  —  1 » i M 1 

1  SUM=aUM*ULl IP i j)*XI J) 

2  X(I)-b  l  i  P)  -0DM 
IP-  IPO  IN) 

X{f\J=X(j\)/GLllP,NI 
LO  4  I  DACK-2  » i\ 

I-NPI-IdAlX 

ip-  ipsm 

ipi-i+i 
SUM-0 . 0 
CO  3  J— IPijN 

3  SGM=oUM+oL{IP»J)*X(J) 

4  xm  =  (xm-3oM)/0LliP,I) 

KtTUKN 

END 
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S  0  d  R  0  u  7  iNt  I  tfPKU  V  (  l\  f  A  ,  UL  ,  d  »  X  ,  0  i  G  1  T  3  7 

Oi  MEi\S  ION  A ( 20 , 20 i  ,  01 { 20 , 20 )  »  6  l  20  )  ,  X  t  2  J  )  ,  K l 20 )  ,  Ox  (  20  J 

OOuriLt  PRECldlUiN  3  u  M  »  X  j  »  A  I  j 

tP  3=  1 .OE-0 

IT.MAX=  12 

X IM  C  K  M  =  C  .  C 

00  i  1=1, l\ 

1  XNGRM=AMAxl(  XMUkM,  AG3UU  )  )  7 
IP  I  Xi\U KjM  73,2,3 

2  OIGI r3=-ALuG10l EPS) 

GO  TU  1C 

3  00  9  I  T tR= 1 ,  I TMAX 
00  5  i  =  1  ,  i\ 

3uM=0. 0 

00  A  J  =  1  ,  i\ 

X J=X I J 7 
A 1  J  =  A  (  i  ,  j  j 

4  3UM=3oRt A I 0*Xj 
30M=  d l  I 7  -SuiM 

5  R  (  1  7  =  3  0  M 

CALL  oLLVtlNfUL|K,uXJ 
u  X  R  C  K  M  =  C  .  C 
00  6  l  =  l,l\ 

T  =  Xl  1  7 

X(  IJ =X  i  I  ) +0 X {  i  7 

OXNuRiv!  =  ARAXi  l  OXNGRiM,  Abo ( X l  i  )  ~T  7  j 
o  CURT  1  NOE 

IF  4  iTtR-17 3, 7,d 

7  ClGiTS  =  -AL0G10lANAXlICXRCRM/XRL,RM,LP3)  7 
3  IF ( OXNlAM-EP 3*XRUkM7  10  ,  lb  ,  0 
O  CORTINCc 

CALL  3  i Ro  l  3  7 
10  RETCRR 
ERG 


' 
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3  0  b  K  u  U  7  i  i\  c  i  l\  V  L  K  T  (N»A,AlNVi 

DIMENSION  A  (  2u  ,  20  i  ,  A  I NV l  20  *  20  )  ,  0  (  20)  ,  Ju l A 0 , 20 )  ,  a  (  20  i 
u  A  L  L  G  t  u  u  ivi  A  i  i\|  A(  U  L  »  U  L  T  >  1  5  GN  U  ) 

Cu  3  J=1,N 
CG  1  1  =  1, N 

b(  n  =  o  .0 

I  )  =  1  •  0 

1  CON  TINGE 

GALL  SGLVt(N,ULib|X) 

GALL  IMPKUtflN, A, UL,b,X, DIGITS! 

Du  2  I  —  1  » i\ 

2  a i N V ( I ,  J)  =  X(  IJ 

3  CUNT i N U E 
RETURN 
END 


b  U  d  R  o  U  T  I  N  e  M  I  N  Pi  A  A  I  Fl » l\  t  A  »  d  »  IV  jX  J 


M  IS  The  N U M ci h K  lF  EQUATIONS. 

N  I  b  I  h  c  NuMd  eR  up  Ui\  K  h  C  w  N  S  • 

A  I  b  The  Pi  -  b  Y  ~  N  MATRIX  CP  COEFFICIENTS# 
d  lb  ThE  KlGHI-hANJ  biCe* 

M  t  N  i  A »  6  MCbT  d  E  buPPLiEb  o  Y  THE  L  A  L  L  i  i\  o  PRUoRAM. 
IV  lb  THe  VECTOR  uF  I  No  I C Ed  uF  Trie  EwUATIUNS 
PREdEhToY  UNGER  OuNd  I  uERm  T 1  Ui\  •  Ah  I N  i  T  1  A  i_ 
o  U  E  d  d  M  0  b  I  bt  dUPPLlLO  d>  ThE  CALLING 
PRCGRaR  •  LN  fseTuKi^iiV  CONTAINS  FHe  iNUlCEb 
CiF  The  EUUATIUNb  uF  The  CRITICAL  dUdbVbiLM. 

X  lb  AN  ( N'T  1 )  ~CCMPGNEnT  VEoTuR.On  KeTuRN» 

x(1)i...x(n)  c  u  n  t  a  i  n  the  Pi  i  n  i  pi  a  x  d  u  l  u  r  I  U  N  » 

An C  X  (  N  + 1 )  C C i\ T  A I  Nb  The  PliNIiTAA  e RK Or\  . 

GlMtNdluN  A  l 100 » i  9  J  ,b(  100)  » I  V ( 2  G  J  »  A  (  2  0  J 
UlMENdiCN  A  A  (  2  0  »  L  0  )  >  u 1 eO  >  20 )  »ULl2CfbO) 

REAL  b  (  20  )  ,  Rhu  ( 100  )  ,  LAiMua  (  CO  ) 

LHJUBlE  PRECIdiuN  Ci\J  »  u  J  i\  j  AlJ  jOIKjXjiLJ  »bUM 

b  T  EP  1 

N  P  1  =  N  + 1 
CO  10  K  =  1 »  N  P  1 
Lu  5  I  —  1  >  i\ 

I F {  I-K) i,e,2 

1  IR=1V(I) 

GC  I C  3 

2  IR=IVlI+l) 

2  CU  A  J  =  i  i  N 

h  AA(I»JI  —  A l Ik i j  J 
d  CuNT  1NCE 

CALL  bEtuMP  (  i\  j  aA  j  UL  »  bt  1  t  IbuNUl 
bTEP  2 

10  SlK^IboNC^I-lJ^^K 
b  I  EP  3 

LG  20  I -=  1 ,  IN! P  1 
I R-  I  V  (  I  ) 

CU  13  J-  1  t N 
Id  AALi»Jj=AlIKjJ) 

20  AAl I i NPl ) =d ( I ) 

call  Invert i n  p 1 »  a  a  »  d ) 


bTEP  A 


o  o  o  o 
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2  5  UU  3  0  K=lfi\Pl 
S  u  ivi  =  0 . 0 
UU  2  9  J=lfiMPl 
IR-  I  V ( o  ) 

UK j=U l K, JJ 
d I R- 3  l  iK) 

29  SUM  =  SUM  +  UKj*diR, 

30  X { K ) =  S  u M 
C 

C  4  STEP  5 

C 

1ALPHA-1 
RhUM A A  =  0  . 0 
b  u  40  i =  1 1 A 
SUM=0.C 
Uu  35  J  —  1  »  im 
A  I 0=  A {  L , Ji 
Au=X  I  J  ) 

3d  SUM=SuM+AIJ$Xj 
r.HUii)=diiJ-^Ui'/l 

d  T  EP  o 

I R  >-*  u  :>  (RHLl  I  )  )  —  K  h  U  f-i  A  a  )  4  0  ?  4  0  »  3  5 
39  RHLRAX=AbS(RhL(iJ ) 

I  ALPHA=  I 
AO  COiM  T  i  IMUt 

STEP  7 

IR  RHUMAX-X  (NPiJ  )  /  1,71,42 
42  Cu  43  K  =  1  » IMP  1 

1 F (  lALPHA-IV  IKJ  J  4  3 >  7  1 »  4  3 
4  3  CbKTIfMLE 

STEP  8 

0=  3UIM  l  Khc  (  i  ALPHA  I  ) 

STEP  9 

RMAX=-1E75 
I  K=  1 

CU  50  K-l»i\Pl 
SUM-0 .  G 
Uu  4  D  J—  I  »  N 
U  J  K  —  U  i  >J  »  K  ) 

4  5  SUM-SUM  +  A I  I  ALPHA i J )  *UJk 
L A  M  U  A ( K ) =SUM  +  U*D<NP1 (KJ 
IF(U(fMPl»KJ  j  4  7  i  5  5  »  4  7 


C 

c 

c 


u 

c 

b 


c 


n  o  o  0-0  0 
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C 

C  6T£P  10 

C 

4  7  T=LAMuAIk)*u/D(NP1»n) 

I  P  l  I  -  k  M  A  X  )  o  C  >  5  0  »  4  9 
4  9  RM AX=  T 
1  R  =  K 

50  Coin T  i NU £ 

1  F  (  L  AMo  A  (  I  r\  )  )59|5Df59 
55  0  A  £  L  a  I NG  1  4  ) 

5  T  t  P  11 

d  9  DO  60  K=l»i\Pl 

oO  D(K»liO  —  J  (  n  »  IKJ/D6LE1LAHJA(  iKi  J 
DO  70  K-=  1 » l\P  1 
CO  7C  J=1»inP1 
!F(J-lKlo9»7Cfo9 

69  OK  j  =  u  ( is  »  J  ) 

LJ=LAMi>AlJl 
ClK»vj)=Oisj~LJ:^L'lKf  IK  J 

70  C  u  l\  T  1  N  0 1 

5T£P  12 

IV  l  UJ  =1  ALPHA 
GO  T c  2  5 

71  R  £  T  U  K  N 
£i\iC 
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b U  b KU U  T  I  IN  E  K  E  M  t  g  l  M  M »  N  t  P T S  »  I  V  »  C  ) 

C 

C  Mi 'l  1 5  The  N  U  M  d  E  K  C  F  Pei  MTS  IN  I  H  E  INITIAL  DIoCkcTc 
C  P  C  I N  T  SET. 

C  N  IS  THl  NUMdEK  CP  AP  PRGX  I  M  AT  In  b  FGnCT  1  GNb  • 

L  P  T  d  IS  TNt  v  t  L  T  l  k  uP  .wodCIbbAb  LP  THE  Pul  Nib  IN 
u  T  Pi  L  uISCklT  t  PCINI  b  c  T  • 

C  Iv  IS  THE  VllIuK  uh  INDICES  GP  PG I  NTS  IN  Thl 
C  CkITICAl  bUdSYSTEM  RETURNED  dY  MlNMAX. 

l  C  IS  An  l  N  +  1 J -CGMPGNt  i\T  v  cC  Tun  •  C  I  I )  »  .  •  •  C  (  N  ) 

C  C  U  NT  A  I N  I H  E  M I N I M  I  L I N  u  CUcFPICIENTb.  u  (  N  +  1 ) 

C  LGKTAINS  The  ASSGCImTEu  MINIMAX  EkKuR. 

C 

DIMENdiuN  P  I  S  1 1 0  G )  tiv(20}#Cl20)»A{lCG»ldjfDlluGJ 
M=MP 

EP  S= i . CE- j 
I TM A  A=  9C 
C 

C  EPb  AND  ITMAa  APE  CHUSEn  uY  Irit  PkUgKAMMcK • 

C  SuT  gP  IHe  M  A  T  k  I  a  a  AND  ThE  VclTGin  b  Purs  M I N  M  A  A  • 

C 

NP  1  =  N+1 
DU  5  I  =  1 1  K 
X=  P  T  b  l  II 
Y  =  bCuf\D  l  X ) 

B(  I )=P ( A,  Y ) 

DC  d  J - i j  N 
d  A(Ijj)-blJfXjY) 

DU  2d  I T EK= 1 » I TMAX 
C 

C  SELECT  Cu  EPF  I C  I  EM  b  T  u  Pi  1 1  ^  I M I  L  E  MAXIMUM  c  k  K  U  k 
D  UN  DIglKlTc  P  u  I N  T  bET  PRESENTLY  dlINg 

C  CCNSI LEKcD. 

c 

call  MINMAXiP  ,Nf A,d,  IVtU) 

c 

C  PiND  FLINT  (CAlLLu  XMAXj  CP  MAXIMUM  EKRGX 
C  IuAlLEu  eMAAJ  DblNu  lUcPPIu  IE  in!  b  UaLLULATEG 

C  P R C M  PKEVlGCb  bTAFEMENT . 

C 

U ALL  ERRMAX(N,u,AMAA,£MAX) 

C 

u  cXIT  RCuTINc. 

c 

IP! EMAa-C ( NP  I ) )d0 , dO , Id 
15  IF  I  (EMAX-C(NPi)  )/i_MAX-EPb)dO»dC»16 
C 

C  INTKUuDCE  POINT  Lp  MAXIMUM  ERkGR  INTO  GlbCKETE 
C  PCInT  SET.  CPD a  T  t  a  AND  D. 

C 


16  N.  =  to  +  1 

P  T  S  (  M )  =  X  K  A  X 
YMAX-BLLMj  IXMAXJ 
B(MJ=MaMmX»  Y  iMaX  ) 

U  G  2  0  J -  1  }  N 

20  A (M  ,  J) =G ( J , XPAX, YMAX J 
2d  c g in r i o t 

CALL  SINGId) 

60  8 1 TU  KN 
tNti 


r  r 
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o  U  d  K  J  U  l  I  iN  b  b  K  K  M  A  X  l  k »  C  »  A  M  A  A  »  b  M  A  A  ) 

C 

u  A  A  i\U  b  ucNuTc  ihL  bi\U  P  □  I N  !  0  UF  The  k  A  k  c  u  F  T  M  l 
C  iNUtPbkUbkr  VhKI  Able.  TFIcV  Pi  U  G  T  be  GJPPuluJ 

u  bY  TFic  pKuuHAMMbk. 

C 

DIMbkSiOk  C  l  2  0  ) 

PI  =3  .  i Al  GS3 
k2C=20*k 
bMAX-0 .0 
Du  10  is—  1»N2C 
C 

U  TFib  K  T  Ft  Z  t  K  U  UF  Thb  u  F  t  u  Y  Sriev'  PULYkuMi.Au  UF 

JbGkub  20^k » kukMAL  1  ZEu  Tu  Trie  IkTcKVAb  (A»bJ» 

X—  C  •  3  *  i  l  b  -  A  }  *  C  U  G  IP  i  *  ( i\  -  0  »b)/k20J  +  d  +  A  i 
Y  =  bCui\u(X) 

E=At3GlriX»Y ) - u l k  »  0  »  X  »  Y  J J 
IFlbivAX-t)b»  IGi  10 
3  b  iM  A  X  =  c 
aMAX^X 
1C  CUM  I N  U  b 
KC  I  UKIM 
bkD 


F u K C T  I  L  i\  U(i'ifLF^»AjY) 

U  I  ivi  1 1\  S  i  u  K  C  F  S  l  2  0  ) 
bUULLt  P  K  c C  I  b l  J in  SUM, OK 
S  u  M  =  0 . 0 
LU  10  A  =  1  ,  N 
Ck=CF  S l K ) 

10  SUM-  SUtt+Cft*U  (  K  ,  X  ,  Y  ) 
u=  Sum 
H 1 1 U  K  M 
ti\U 


FUNCTION  bGNlZ) 
IF  l  Z  )  i  »  2  »  3 

1  S  u  N  -  -  1 .  b 
RETURN 

2  SoN=C.C 
h  t  T  U  kN 

3  SbN=  1  •  C 
return 
end 
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FUNCTICN  FIX»YJ 

F-ALUb ( t  X  >  1  . b ) v*2  +  i 

Rt I  URN 

END 


FUi\CT  I  UN  G  t  in  »  X  #  Y  ) 
COMPLEX  L 
KM  1  =  ft"  i 
IP  t  K M  1  }  1»  1 »  2 

1  G= 1  •  0 
RETCKN 

2  Z=CMPLX  (  X  »  Y  ) 
G=kEAL 1 ) 
RETCKN 

EMU 
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f  Ut\C  <  I  Ci\  tiCUiNL  1  X  ) 
dbUNu-C  •  i  1  ~  X  *  ^  2  ) 

H 1 1 U  H  N 
cNU 


SboRLUT 1 U  £  iiNollnHY) 

U  u  T  U  {l»Ef3j4»3)fiWFiY 
1  hKI TE { b i i  1 ) 

R  E  T  G  K  i\ 


2  aR17c(g,U) 
h  E  I  L  K  l\ 

3  rtRITtlfc, 13) 
K  £  T  L.  Kl\ 


4  «KITEto,  i^J 
K  E  T  0  R  i\ 


3  ftKlTc{6,13) 
k£TL«l\ 


11  FlkMATI'C 
1  £  F  u  K  H  A  T  (  '  0 

13  FURR AT  l  1 C 

14  FURR AT ( ' C 
1 3  FuRMaT l 1 0 

ENl. 


A  T  K  i  X  vilTh  *l£KU  R  Q  irt  I  i\  LitCOMP.  1  ) 


n 

ii  iWulak  MmT^IX  In  JcCUMP  •  '  ) 
l\G  CuN V £ Ku ci'io L  in  IRiFkUV.’J 
HrtAK  UUimU  I  T  i  Ul\l 
i^U  C  UU  V  £  RU  tl'Uc 


v'IUlaTEl)  1  i\ 
I  U  KcMLo.  ‘  ) 


M  i  l\  i'l  A  A  • 
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Minimax  Coefficients  for  the  Bean  Problem 


r 


(x,y)  =  < 

0 ,  x  <  0 

x,  x  >  0 

n  =  3  : 

.132231 

.  461538 

LT\ 

II 

£ 

.  093^04 

.486649 

n  =  7  : 

.092597 

.536986 

. 180471 

. 441858 

n  =  19 : 

.075025 

. 514878 

-.224194 

1.03329 

-2 . 60345 

6.39725 

-9.69454 

-.596876 

.092308 
. 008827 

.608221 

-.034363 

-.105155 

.653991 

.018744 

-.209529 

.989064 

.015350 

-2 .27918 

-.035892 

. 213942 

-.026738 

1.27328 

4.87422 

11.2224 

-5 . 88364 

f(x,y)  =  < 

r 

0  ,  x  <  0 

x3  x  >  0 

n  =  3  : 

.031894 

. 120151 

-.000872 

n  =  5  : 

.017586 

. 014766 

-.014176 

. 199221 

- .  006899 

n  =  7  : 

.013713 

.089248 

-.032558 

. 196624 

.  034272 

.  212272 

.  065763 

1 — 1 

II 

£ 

.010926 

.079131 

-.048096 

.236794 

-.024109 

.241918 

. 200170 

-.175985 

.255733 

-.415811 

-.305695 

.488001 

-.453883 

1.14545 

.  212120 

-1.02356 

.  561388 

-1.88455 

.951538 
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f(x,y ) 


x  <  0 
x  >  0 


n  =  3  : 

.  007117 

.032789 

n  =  5  : 

.  003361 

.  023898 

n  =  7  : 

.002812 

.021016 

.055872 

.019046 

n  =  19 : 

.  002201 

.016757 

.072625 

.045666 

-.035986 

.078506 

-.229025 

.153825 

0.0 

-.002281 

.  056748 

. 005689 

-.007908 

.049424 

.011138 

-.009721 

.055003 

- .  004772 

. 002280 

.071505 

- . 066786 

-.096603 

. 184970 

.009972 

-.3^7786 

.244229 

<< 

11 

X 

e  cos  y  + 

log[ ( 1-y )2 

+  x2] 

n  =  3  : 

1.19201 

.  861679 

-1.62099 

n  =  5  : 

1.04665 

.993063 

-1. 80887 

1.40812 

-.222934 

n  =  7  : 

1.01423 

.  991961 

-1.94136 

1.36212 

-.017493 

.  309424 

. 610016 

n  =  19 : 

1.00007 

.999835 

-1.99968 

1.49914 

-.000545 

. 168320 

.664599 

-.453063 

.  004362 

-.000485 

-.388996 

.309699 

-.018050 

.039101 

.243903 

-.168254 

.007315 

-.114652 

- . 120182 

. 

